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:
James E McClure 2021-03-31 13:03:56 -04:00
commit dd94f96215
11 changed files with 8748 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
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 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); 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 // ION TRANSPORT MODEL
extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *dist, double *Den, int start, int finish, int Np); extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *dist, double *Den, int start, int finish, int Np);
@ -193,11 +206,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, 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); 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, 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); 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, 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); double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np);
@ -214,6 +228,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, 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); 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, 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); 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

1304
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

1340
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; Fx = Fy = Fz = 0.0;
gamma=1e-3;//surface tension gamma=1e-3;//surface tension
W=5.0;//interfacial thickness W=5.0;//interfacial thickness
beta = 12.0*gamma/W; //beta = 12.0*gamma/W;
kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma //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; Restart=false;
din=dout=1.0; din=dout=1.0;
flux=0.0; flux=0.0;
@ -136,8 +138,10 @@ void ScaLBL_FreeLeeModel::ReadParams(string filename){
outletA=0.f; outletA=0.f;
outletB=1.f; outletB=1.f;
//update secondary parameters //update secondary parameters
beta = 12.0*gamma/W; //beta = 12.0*gamma/W;
kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma //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) //if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details)
BoundaryCondition = 0; BoundaryCondition = 0;
@ -795,27 +799,31 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
// Compute the Phase indicator field // Compute the Phase indicator field
// Read for hq happens in this routine (requires communication) // Read for hq happens in this routine (requires communication)
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMAL 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->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); 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 // Perform the collision operation
// Halo exchange for phase field //ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastInterior(), Np); //ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_Comm_WideHalo->Send(Phi); //ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm_WideHalo->Recv(Phi); ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL
if (BoundaryCondition > 0 && BoundaryCondition < 5){ if (BoundaryCondition > 0 && BoundaryCondition < 5){
//TODO to be revised //TODO to be revised
// Need to add BC for hq!!! // Need to add BC for hq!!!
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
} }
// Halo exchange for phase field
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL ScaLBL_Comm_WideHalo->Send(Phi);
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, 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); 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->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
// Set BCs // Set BCs
@ -832,7 +840,9 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(gqbar); 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); kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
@ -841,24 +851,29 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
timestep++; timestep++;
// Compute the Phase indicator field // Compute the Phase indicator field
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMA 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->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); 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 // 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_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); //ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_Comm_WideHalo->Send(Phi); //ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm_WideHalo->Recv(Phi); ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL
if (BoundaryCondition > 0 && BoundaryCondition < 5){ if (BoundaryCondition > 0 && BoundaryCondition < 5){
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
} }
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL // 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, //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); 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->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
// Set boundary conditions // 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_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); kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
//************************************************************************ //************************************************************************
@ -1055,6 +1072,13 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){
fwrite(PhaseField.data(),8,N,PFILE); fwrite(PhaseField.data(),8,N,PFILE);
fclose(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); ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
FILE *VELX_FILE; FILE *VELX_FILE;
sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank); 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): 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), 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) 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; REVERSE_FLOW_DIRECTION = false;
@ -43,6 +43,8 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
Restart=false; Restart=false;
din=dout=1.0; din=dout=1.0;
flux=0.0; flux=0.0;
RecoloringOff = false;
W=1.0;
// Color Model parameters // Color Model parameters
if (greyscaleColor_db->keyExists( "timestepMax" )){ if (greyscaleColor_db->keyExists( "timestepMax" )){
@ -85,6 +87,12 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
if (greyscaleColor_db->keyExists( "flux" )){ if (greyscaleColor_db->keyExists( "flux" )){
flux = greyscaleColor_db->getScalar<double>( "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; inletA=1.f;
inletB=0.f; inletB=0.f;
outletA=0.f; outletA=0.f;
@ -212,6 +220,7 @@ void ScaLBL_GreyscaleColorModel::ReadInput(){
} }
void ScaLBL_GreyscaleColorModel::AssignComponentLabels() void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
{ {
// Initialize impermeability solid nodes and grey nodes // 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]); //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){ if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx]; AFFINITY=AffinityList[idx];
label_count[idx] += 1.0; label_count[idx] += 1.0;
idx = NLABELS; idx = NLABELS;
//Mask->id[n] = 0; // set mask to zero since this is an immobile component //Mask->id[n] = 0; // set mask to zero since this is an immobile component
@ -575,6 +585,71 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
delete [] Permeability; delete [] Permeability;
} }
void ScaLBL_GreyscaleColorModel::AssignGreyscalePotential()
{
double *psi;//greyscale potential
psi = new double[N];
size_t NLABELS=0;
signed char VALUE=0;
double AFFINITY=0.f;
auto LabelList = greyscaleColor_db->getVector<int>( "ComponentLabels" );
auto AffinityList = greyscaleColor_db->getVector<double>( "ComponentAffinity" );
NLABELS=LabelList.size();
//first, copy over normal phase field
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
VALUE=id[n];
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
idx = NLABELS;
}
}
// fluid labels are reserved
if (VALUE == 1) AFFINITY=1.0;
else if (VALUE == 2) AFFINITY=-1.0;
psi[n] = AFFINITY;
}
}
}
//second, scale the phase field for grey nodes
double Cap_Penalty=1.f;
auto GreyLabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
auto PermeabilityList = greyscaleColor_db->getVector<double>( "PermeabilityList" );
NLABELS=GreyLabelList.size();
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
VALUE=id[n];
Cap_Penalty=1.f;
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
if (VALUE == GreyLabelList[idx]){
Cap_Penalty=alpha*W/sqrt(PermeabilityList[idx]/Dm->voxel_length/Dm->voxel_length);
idx = NLABELS;
}
}
//update greyscale potential
psi[n] = psi[n]*Cap_Penalty;
}
}
}
ScaLBL_CopyToDevice(Psi, psi, N*sizeof(double));
ScaLBL_Comm->Barrier();
delete [] psi;
}
void ScaLBL_GreyscaleColorModel::Create(){ void ScaLBL_GreyscaleColorModel::Create(){
/* /*
* This function creates the variables needed to run a LBM * This function creates the variables needed to run a LBM
@ -595,6 +670,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
// ScaLBL_Communicator ScaLBL_Comm(Mask); // original // ScaLBL_Communicator ScaLBL_Comm(Mask); // original
ScaLBL_Comm = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask)); 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 = 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; int Npad=(Np/16 + 2)*16;
if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); 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 **) &Bq, 7*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size); ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nx*Ny*Nz); 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 **) &Pressure, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
//ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np); //ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np);
@ -667,6 +744,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
AssignComponentLabels();//do open/black/grey nodes initialization AssignComponentLabels();//do open/black/grey nodes initialization
AssignGreySolidLabels(); AssignGreySolidLabels();
AssignGreyPoroPermLabels(); AssignGreyPoroPermLabels();
AssignGreyscalePotential();
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity); Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity);
ScaLBL_Comm->RegularLayout(Map,Porosity_dvc,Averages->Porosity);//porosity doesn't change over time ScaLBL_Comm->RegularLayout(Map,Porosity_dvc,Averages->Porosity);//porosity doesn't change over time
} }
@ -931,9 +1009,11 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Read for Aq, Bq happens in this routine (requires communication) // Read for Aq, Bq happens in this routine (requires communication)
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL 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_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->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); 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 // Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
@ -943,14 +1023,20 @@ void ScaLBL_GreyscaleColorModel::Run(){
} }
// Halo exchange for phase field // Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi); ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4 ScaLBL_Comm_Regular_2->SendHalo(Psi);
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, 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 ////Model-2&3
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity, //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, // 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, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular_2->RecvHalo(Psi);
ScaLBL_Comm_Regular->RecvHalo(Phi); ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
@ -968,10 +1054,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
} }
//Model-1&4 //Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure, 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, 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 ////Model-2&3
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity, //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, // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
@ -983,9 +1073,11 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Compute the Phase indicator field // Compute the Phase indicator field
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL 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_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->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); 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 // Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL 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->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
} }
ScaLBL_Comm_Regular->SendHalo(Phi); ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4 ScaLBL_Comm_Regular_2->SendHalo(Psi);
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, 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 ////Model-2&3
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity, //ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, // 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, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular_2->RecvHalo(Psi);
ScaLBL_Comm_Regular->RecvHalo(Phi); ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier(); ScaLBL_Comm->Barrier();
@ -1020,10 +1118,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
} }
//Model-1&4 //Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure, 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, 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 ////Model-2&3
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity, //ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, // rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
@ -1575,6 +1677,13 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
fwrite(PhaseField.data(),8,N,OUTFILE); fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(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); ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
FILE *AFILE; FILE *AFILE;
sprintf(LocalRankFilename,"A.%05i.raw",rank); sprintf(LocalRankFilename,"A.%05i.raw",rank);

7
models/GreyscaleColorModel.h
View File

@ -39,6 +39,8 @@ public:
double Fx,Fy,Fz,flux; double Fx,Fy,Fz,flux;
double din,dout,inletA,inletB,outletA,outletB; double din,dout,inletA,inletB,outletA,outletB;
double GreyPorosity; 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 Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs; int rank,nprocx,nprocy,nprocz,nprocs;
@ -48,7 +50,8 @@ public:
std::shared_ptr<Domain> Mask; // this domain is for lbm std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm; std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular; std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
std::shared_ptr<GreyPhaseAnalysis> Averages; std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular_2;
std::shared_ptr<GreyPhaseAnalysis> Averages;
// input database // input database
std::shared_ptr<Database> db; std::shared_ptr<Database> db;
@ -70,6 +73,7 @@ public:
double *Pressure; double *Pressure;
double *Porosity_dvc; double *Porosity_dvc;
double *Permeability_dvc; double *Permeability_dvc;
double *Psi;
private: private:
Utilities::MPI comm; Utilities::MPI comm;
@ -86,6 +90,7 @@ private:
void AssignComponentLabels(); void AssignComponentLabels();
void AssignGreySolidLabels(); void AssignGreySolidLabels();
void AssignGreyPoroPermLabels(); void AssignGreyPoroPermLabels();
void AssignGreyscalePotential();
void ImageInit(std::string filename); void ImageInit(std::string filename);
double MorphInit(const double beta, const double morph_delta); double MorphInit(const double beta, const double morph_delta);
double SeedPhaseField(const double seed_water_in_oil); double SeedPhaseField(const double seed_water_in_oil);

2
tests/lbpm_freelee_simulator.cpp
View File

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