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

Merge branch 'master' into greyscale_dev

Browse Source
This commit is contained in:
JamesEMcclure 2021-06-16 13:45:13 -04:00
commit b4f4607db0
46 changed files with 1561 additions and 678 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists.txt
View File

@ -18,8 +18,7 @@ SET( TEST_MAX_PROCS 16 )
# Initialize the project
PROJECT( ${PROJ} LANGUAGES CXX )
PROJECT( ${PROJ} LANGUAGES CXX)
# Prevent users from building in place
IF ("${CMAKE_CURRENT_SOURCE_DIR}" STREQUAL "${CMAKE_CURRENT_BINARY_DIR}" )

49
analysis/ElectroChemistry.cpp
View File

@ -53,19 +53,32 @@ void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poiss
Poisson.getElectricPotential(ElectricalPotential);
/* local sub-domain averages */
double rho_avg_local[Ion.number_ion_species];
double rho_mu_avg_local[Ion.number_ion_species];
double rho_mu_fluctuation_local[Ion.number_ion_species];
double rho_psi_avg_local[Ion.number_ion_species];
double rho_psi_fluctuation_local[Ion.number_ion_species];
double *rho_avg_local;
double *rho_mu_avg_local;
double *rho_mu_fluctuation_local;
double *rho_psi_avg_local;
double *rho_psi_fluctuation_local;
/* global averages */
double rho_avg_global[Ion.number_ion_species];
double rho_mu_avg_global[Ion.number_ion_species];
double rho_mu_fluctuation_global[Ion.number_ion_species];
double rho_psi_avg_global[Ion.number_ion_species];
double rho_psi_fluctuation_global[Ion.number_ion_species];
double *rho_avg_global;
double *rho_mu_avg_global;
double *rho_mu_fluctuation_global;
double *rho_psi_avg_global;
double *rho_psi_fluctuation_global;
for (int ion=0; ion<Ion.number_ion_species; ion++){
/* local sub-domain averages */
rho_avg_local = new double [Ion.number_ion_species];
rho_mu_avg_local = new double [Ion.number_ion_species];
rho_mu_fluctuation_local = new double [Ion.number_ion_species];
rho_psi_avg_local = new double [Ion.number_ion_species];
rho_psi_fluctuation_local = new double [Ion.number_ion_species];
/* global averages */
rho_avg_global = new double [Ion.number_ion_species];
rho_mu_avg_global = new double [Ion.number_ion_species];
rho_mu_fluctuation_global = new double [Ion.number_ion_species];
rho_psi_avg_global = new double [Ion.number_ion_species];
rho_psi_fluctuation_global = new double [Ion.number_ion_species];
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
rho_avg_local[ion] = 0.0;
rho_mu_avg_local[ion] = 0.0;
rho_psi_avg_local[ion] = 0.0;
@ -90,7 +103,7 @@ void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poiss
}
}
for (int ion=0; ion<Ion.number_ion_species; ion++){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
rho_mu_fluctuation_local[ion] = 0.0;
rho_psi_fluctuation_local[ion] = 0.0;
/* Compute averages for each ion */
@ -108,7 +121,7 @@ void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poiss
if (Dm->rank()==0){
fprintf(TIMELOG,"%i ",timestep);
for (int ion=0; ion<Ion.number_ion_species; ion++){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
fprintf(TIMELOG,"%.8g ",rho_avg_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_mu_avg_global[ion]);
fprintf(TIMELOG,"%.8g ",rho_psi_avg_global[ion]);
@ -147,7 +160,7 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
visData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
auto ElectricPotential = std::make_shared<IO::Variable>();
std::vector<shared_ptr<IO::Variable>> IonConcentration;
for (int ion=0; ion<Ion.number_ion_species; ion++){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
IonConcentration.push_back(std::make_shared<IO::Variable>());
}
auto VxVar = std::make_shared<IO::Variable>();
@ -163,8 +176,8 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
}
if (vis_db->getWithDefault<bool>( "save_concentration", true )){
for (int ion=0; ion<Ion.number_ion_species; ion++){
sprintf(VisName,"IonConcentration_%i",ion+1);
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
sprintf(VisName,"IonConcentration_%zu",ion+1);
IonConcentration[ion]->name = VisName;
IonConcentration[ion]->type = IO::VariableType::VolumeVariable;
IonConcentration[ion]->dim = 1;
@ -199,8 +212,8 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
}
if (vis_db->getWithDefault<bool>( "save_concentration", true )){
for (int ion=0; ion<Ion.number_ion_species; ion++){
sprintf(VisName,"IonConcentration_%i",ion+1);
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
sprintf(VisName,"IonConcentration_%zu",ion+1);
IonConcentration[ion]->name = VisName;
ASSERT(visData[0].vars[1+ion]->name==VisName);
Array<double>& IonConcentrationData = visData[0].vars[1+ion]->data;

3
analysis/FreeEnergy.cpp
View File

@ -47,8 +47,6 @@ void FreeEnergyAnalyzer::SetParams(){
void FreeEnergyAnalyzer::Basic(ScaLBL_FreeLeeModel &LeeModel, int timestep){
int i,j,k;
if (Dm->rank()==0){
fprintf(TIMELOG,"%i ",timestep);
/*for (int ion=0; ion<Ion.number_ion_species; ion++){
@ -78,7 +76,6 @@ void FreeEnergyAnalyzer::Basic(ScaLBL_FreeLeeModel &LeeModel, int timestep){
void FreeEnergyAnalyzer::WriteVis( ScaLBL_FreeLeeModel &LeeModel, std::shared_ptr<Database> input_db, int timestep){
auto vis_db = input_db->getDatabase( "Visualization" );
char VisName[40];
std::vector<IO::MeshDataStruct> visData;
fillHalo<double> fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1);

20
analysis/morphology.cpp
View File

@ -120,6 +120,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
sendtag = recvtag = 7;
int ii,jj,kk;
int imin,jmin,kmin,imax,jmax,kmax;
int Nx = nx;
int Ny = ny;
int Nz = nz;
@ -128,17 +129,13 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
double void_fraction_new=1.0;
double void_fraction_diff_old = 1.0;
double void_fraction_diff_new = 1.0;
// Increase the critical radius until the target saturation is met
double deltaR=0.05; // amount to change the radius in voxel units
double Rcrit_old;
int imin,jmin,kmin,imax,jmax,kmax;
if (ErodeLabel == 1){
VoidFraction = 1.0 - VoidFraction;
}
// Increase the critical radius until the target saturation is met
double deltaR=0.05; // amount to change the radius in voxel units
double Rcrit_old = maxdistGlobal;
double Rcrit_new = maxdistGlobal;
while (void_fraction_new > VoidFraction)
@ -406,6 +403,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
sendtag = recvtag = 7;
*/
int ii,jj,kk;
int imin,jmin,kmin,imax,jmax,kmax;
int Nx = nx;
int Ny = ny;
int Nz = nz;
@ -417,10 +415,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
// Increase the critical radius until the target saturation is met
double deltaR=0.05; // amount to change the radius in voxel units
double Rcrit_old;
int imin,jmin,kmin,imax,jmax,kmax;
double Rcrit_old = maxdistGlobal;
double Rcrit_new = maxdistGlobal;
//if (argc>2){
// Rcrit_new = strtod(argv[2],NULL);
@ -457,7 +452,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
for (kk=kmin; kk<kmax; kk++){
for (jj=jmin; jj<jmax; jj++){
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(ii,jj,kk) == 2 && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
LocalNumber+=1.0;
@ -578,7 +572,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
// nwp
phase(i,j,k) = -1.0;
}
else{
else{i
// treat solid as WP since films can connect
phase(i,j,k) = 1.0;
}

1
analysis/runAnalysis.cpp
View File

@ -729,7 +729,6 @@ runAnalysis::runAnalysis( ScaLBL_ColorModel &ColorModel)
d_comm = ColorModel.Dm->Comm.dup();
d_Np = ColorModel.Np;
bool Regular = false;
auto input_db = ColorModel.db;
auto db = input_db->getDatabase( "Analysis" );

135
common/ScaLBL.cpp
View File

@ -1046,29 +1046,28 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
}
}
}
int *bb_dist_tmp = new int [local_count];
int *bb_interactions_tmp = new int [local_count];
ScaLBL_AllocateDeviceMemory((void **) &bb_dist, sizeof(int)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &bb_interactions, sizeof(int)*local_count);
int *fluid_boundary_tmp;
double *lattice_weight_tmp;
float *lattice_cx_tmp;
float *lattice_cy_tmp;
float *lattice_cz_tmp;
if(SlippingVelBC==true){
fluid_boundary_tmp = new int [local_count];
lattice_weight_tmp = new double [local_count];
lattice_cx_tmp = new float [local_count];
lattice_cy_tmp = new float [local_count];
lattice_cz_tmp = new float [local_count];
ScaLBL_AllocateDeviceMemory((void **) &fluid_boundary, sizeof(int)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_weight, sizeof(double)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_cx, sizeof(float)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_cy, sizeof(float)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_cz, sizeof(float)*local_count);
}
int *fluid_boundary_tmp;
double *lattice_weight_tmp;
float *lattice_cx_tmp;
float *lattice_cy_tmp;
float *lattice_cz_tmp;
/* allocate memory for bounce-back sites */
fluid_boundary_tmp = new int [local_count];
lattice_weight_tmp = new double [local_count];
lattice_cx_tmp = new float [local_count];
lattice_cy_tmp = new float [local_count];
lattice_cz_tmp = new float [local_count];
ScaLBL_AllocateDeviceMemory((void **) &fluid_boundary, sizeof(int)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_weight, sizeof(double)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_cx, sizeof(float)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_cy, sizeof(float)*local_count);
ScaLBL_AllocateDeviceMemory((void **) &lattice_cz, sizeof(float)*local_count);
local_count=0;
for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){
@ -1081,78 +1080,78 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
neighbor=Map(i-1,j,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i-1) + (j)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 0.0;
}
//}
bb_dist_tmp[local_count++]=idx + 2*Np;
}
neighbor=Map(i+1,j,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i+1) + (j)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 0.0;
}
//}
bb_dist_tmp[local_count++] = idx + 1*Np;
}
neighbor=Map(i,j-1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j-1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = 0.0;
}
//}
bb_dist_tmp[local_count++]=idx + 4*Np;
}
neighbor=Map(i,j+1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j+1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = 0.0;
}
//}
bb_dist_tmp[local_count++]=idx + 3*Np;
}
neighbor=Map(i,j,k-1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = -1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 6*Np;
}
neighbor=Map(i,j,k+1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/18.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 5*Np;
}
}
@ -1170,156 +1169,156 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
neighbor=Map(i-1,j-1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i-1) + (j-1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = 0.0;
}
//}
bb_dist_tmp[local_count++]=idx + 8*Np;
}
neighbor=Map(i+1,j+1,k);
if (neighbor==-1) {
bb_interactions_tmp[local_count] = (i+1) + (j+1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = 0.0;
}
//}
bb_dist_tmp[local_count++]=idx + 7*Np;
}
neighbor=Map(i-1,j+1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i-1) + (j+1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = 0.0;
}
//}
bb_dist_tmp[local_count++]=idx + 10*Np;
}
neighbor=Map(i+1,j-1,k);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i+1) + (j-1)*Nx + (k)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = 0.0;
}
//}
bb_dist_tmp[local_count++]=idx + 9*Np;
}
neighbor=Map(i-1,j,k-1);
if (neighbor==-1) {
bb_interactions_tmp[local_count] = (i-1) + (j)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = -1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 12*Np;
}
neighbor=Map(i+1,j,k+1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i+1) + (j)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 11*Np;
}
neighbor=Map(i-1,j,k+1);
if (neighbor==-1) {
bb_interactions_tmp[local_count] = (i-1) + (j)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = -1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = 1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 14*Np;
}
neighbor=Map(i+1,j,k-1);
if (neighbor==-1) {
bb_interactions_tmp[local_count] = (i+1) + (j)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 1.0;
lattice_cy_tmp[local_count] = 0.0;
lattice_cz_tmp[local_count] = -1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 13*Np;
}
neighbor=Map(i,j-1,k-1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j-1)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = -1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 16*Np;
}
neighbor=Map(i,j+1,k+1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j+1)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = 1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 15*Np;
}
neighbor=Map(i,j-1,k+1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j-1)*Nx + (k+1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = -1.0;
lattice_cz_tmp[local_count] = 1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 18*Np;
}
neighbor=Map(i,j+1,k-1);
if (neighbor==-1){
bb_interactions_tmp[local_count] = (i) + (j+1)*Nx + (k-1)*Nx*Ny;
if(SlippingVelBC==true){
//if(SlippingVelBC==true){
fluid_boundary_tmp[local_count] = idx;
lattice_weight_tmp[local_count] = 1.0/36.0;
lattice_cx_tmp[local_count] = 0.0;
lattice_cy_tmp[local_count] = 1.0;
lattice_cz_tmp[local_count] = -1.0;
}
//}
bb_dist_tmp[local_count++]=idx + 17*Np;
}
}
@ -1329,24 +1328,20 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in
n_bb_d3q19 = local_count; // this gives the d3q19 distributions not part of d3q7 model
ScaLBL_CopyToDevice(bb_dist, bb_dist_tmp, local_count*sizeof(int));
ScaLBL_CopyToDevice(bb_interactions, bb_interactions_tmp, local_count*sizeof(int));
if(SlippingVelBC==true){
ScaLBL_CopyToDevice(fluid_boundary, fluid_boundary_tmp, local_count*sizeof(int));
ScaLBL_CopyToDevice(lattice_weight, lattice_weight_tmp, local_count*sizeof(double));
ScaLBL_CopyToDevice(lattice_cx, lattice_cx_tmp, local_count*sizeof(float));
ScaLBL_CopyToDevice(lattice_cy, lattice_cy_tmp, local_count*sizeof(float));
ScaLBL_CopyToDevice(lattice_cz, lattice_cz_tmp, local_count*sizeof(float));
}
ScaLBL_CopyToDevice(fluid_boundary, fluid_boundary_tmp, local_count*sizeof(int));
ScaLBL_CopyToDevice(lattice_weight, lattice_weight_tmp, local_count*sizeof(double));
ScaLBL_CopyToDevice(lattice_cx, lattice_cx_tmp, local_count*sizeof(float));
ScaLBL_CopyToDevice(lattice_cy, lattice_cy_tmp, local_count*sizeof(float));
ScaLBL_CopyToDevice(lattice_cz, lattice_cz_tmp, local_count*sizeof(float));
ScaLBL_DeviceBarrier();
delete [] bb_dist_tmp;
delete [] bb_interactions_tmp;
if(SlippingVelBC==true){
delete [] fluid_boundary_tmp;
delete [] lattice_weight_tmp;
delete [] lattice_cx_tmp;
delete [] lattice_cy_tmp;
delete [] lattice_cz_tmp;
}
delete [] fluid_boundary_tmp;
delete [] lattice_weight_tmp;
delete [] lattice_cx_tmp;
delete [] lattice_cy_tmp;
delete [] lattice_cz_tmp;
}
void ScaLBL_Communicator::SolidDirichletD3Q7(double *fq, double *BoundaryValue){

5
cpu/BGK.cpp
View File

@ -1,5 +1,4 @@
extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){
int n;
// conserved momemnts
double rho,ux,uy,uz,uu;
// non-conserved moments
@ -111,14 +110,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int
}
extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){
int n;
// conserved momemnts
double rho,ux,uy,uz,uu;
// non-conserved moments
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
int nr1,nr2,nr3,nr4,nr5,nr6,nr7,nr8,nr9,nr10,nr11,nr12,nr13,nr14,nr15,nr16,nr17,nr18;
int nread;
for (int n=start; n<finish; n++){
// q=0
@ -275,4 +272,4 @@ extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int star
rlx*0.02777777777777778*(rho - 3.0*(uy-uz) + 4.5*(uy-uz)*(uy-uz) - uu) - 0.08333333333*(Fy-Fz);
}
}
}

17
cpu/Color.cpp
View File

@ -919,17 +919,14 @@ extern "C" void ScaLBL_D3Q19_ColorCollide( char *ID, double *disteven, double *d
extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A_odd, double *B_even, double *B_odd,
double *Den, double *Phi, double *ColorGrad, double *Velocity, double beta, int N, bool pBC)
{
int n;
char id;
int idx,n,q,Cqx,Cqy,Cqz;
// int sendLoc;
double f0,f1,f2,f3,f4,f5,f6;
double na,nb,nab; // density values
double ux,uy,uz; // flow velocity
double nx,ny,nz,C; // color gradient components
double a1,a2,b1,b2;
double sp,delta;
double delta;
//double feq[6]; // equilibrium distributions
// Set of Discrete velocities for the D3Q19 Model
//int D3Q7[3][3]={{1,0,0},{0,1,0},{0,0,1}};
@ -1255,7 +1252,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, do
double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
int ijk,nn,n;
int ijk,nn;
double fq;
// conserved momemnts
double rho,jx,jy,jz;
@ -1838,7 +1835,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *di
double *Phi, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
int n,nn,ijk,nread;
int nn,ijk,nread;
int nr1,nr2,nr3,nr4,nr5,nr6;
int nr7,nr8,nr9,nr10;
int nr11,nr12,nr13,nr14;
@ -2498,7 +2495,6 @@ extern "C" void ScaLBL_D3Q7_AAodd_Color(int *neighborList, int *Map, double *Aq,
double a1,b1,a2,b2,nAB,delta;
double C,nx,ny,nz; //color gradient magnitude and direction
double ux,uy,uz;
double phi;
// Instantiate mass transport distributions
// Stationary value - distribution 0
for (int n=start; n<finish; n++){
@ -2531,7 +2527,6 @@ extern "C" void ScaLBL_D3Q7_AAodd_Color(int *neighborList, int *Map, double *Aq,
nB = Den[Np + n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
nAB = 1.0/(nA+nB);
Aq[n] = 0.3333333333333333*nA;
Bq[n] = 0.3333333333333333*nB;
@ -2682,7 +2677,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Color(int *Map, double *Aq, double *Bq, doubl
extern "C" void ScaLBL_D3Q7_AAodd_PhaseField(int *neighborList, int *Map, double *Aq, double *Bq,
double *Den, double *Phi, int start, int finish, int Np){
int idx,nread;
int idx,nread;
double fq,nA,nB;
for (int n=start; n<finish; n++){
@ -2842,7 +2837,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_PhaseField(int *Map, double *Aq, double *Bq,
}
extern "C" void ScaLBL_D3Q19_Gradient(int *Map, double *phi, double *ColorGrad, int start, int finish, int Np, int Nx, int Ny, int Nz){
int idx,n,N,i,j,k,nn;
int idx,n,i,j,k,nn;
// distributions
double f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;

20
docs/Makefile Normal file
View File

@ -0,0 +1,20 @@
# Minimal makefile for Sphinx documentation
#
# You can set these variables from the command line, and also
# from the environment for the first two.
SPHINXOPTS ?=
SPHINXBUILD ?= sphinx-build
SOURCEDIR = source
BUILDDIR = $(HOME)/local/doc/build
# Put it first so that "make" without argument is like "make help".
help:
@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
.PHONY: help Makefile
# Catch-all target: route all unknown targets to Sphinx using the new
# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
%: Makefile
@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)

19
docs/README.md Normal file
View File

@ -0,0 +1,19 @@
Dependencies for LBPM documentation
# install sphinx
pip install Sphinx
# foamatting requires sphinx read-the-docs-theme
pip install sphinx-rtd-theme
# equation rendering requires latex and dvipng command
sudo apt-get install dvipng
sudo apt-get install texlive texstudio
sudo apt-get install texlive-latex-recommended texlive-pictures texlive-latex-extra
# To build the docs
Step 1) install dependencies listed above
Step 2) type 'make html' from the docs/ directory
Step 3) point your browser at ~/local/doc/build/html/index.html
#

70
docs/source/conf.py Normal file
View File

@ -0,0 +1,70 @@
# Configuration file for the Sphinx documentation builder.
#
# This file only contains a selection of the most common options. For a full
# list see the documentation:
# https://www.sphinx-doc.org/en/master/usage/configuration.html
# -- Path setup --------------------------------------------------------------
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#
import os
import sys
# sys.path.insert(0, os.path.abspath('.'))
# -- Project information -----------------------------------------------------
project = 'LBPM'
copyright = '2021, James E McClure'
author = 'James E McClure'
# The full version, including alpha/beta/rc tags
release = '1.0'
# -- General configuration ---------------------------------------------------
# Add any Sphinx extension module names here, as strings. They can be
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones.
extensions = [
'sphinx.ext.imgmath'
]
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
# This pattern also affects html_static_path and html_extra_path.
exclude_patterns = []
# -- Options for HTML output -------------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
#
html_theme = 'alabaster'
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
## Read the docs style:
if os.environ.get('READTHEDOCS') != 'True':
try:
import sphinx_rtd_theme
except ImportError:
pass # assume we have sphinx >= 1.3
else:
html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
html_theme = 'sphinx_rtd_theme'
#def setup(app):
# app.add_stylesheet("fix_rtd.css")

20
docs/source/developerGuide/buildingModels/overview.rst Normal file
View File

@ -0,0 +1,20 @@
===========================
Implementing a new LB model
===========================
While LBPM includes a range of fully-functioning lattice Boltzmann models, the commonly used
Bhatnager-Gross-Krook (BGK) model has been deliberately excluded. While the physical limitations
of this model are well-known, implementing the BGK model is an excellent way to understand
how to implement new LB models within the more general framework of LBPM. In this excercise
you will
* learn "what goes where"
* don't modify core data structures (unless you have a really good reason)
* re-use existing components whenever possible

6
docs/source/developerGuide/designOverview/dataStructures.rst Normal file
View File

@ -0,0 +1,6 @@
===============
Data Structures
===============
LBPM includes a variety of generalized data structures to facilitate the implementation
of different lattice Boltzmann models.

18
docs/source/developerGuide/index.rst Normal file
View File

@ -0,0 +1,18 @@
###############################################################################
Developer guide
###############################################################################
The LBPM developer guide provides essential information on how to add new physics
into the framework.
.. toctree::
:glob:
:maxdepth: 2
designOverview/*
buildingModels/*
testingModels/*

9
docs/source/developerGuide/testingModels/unitTests.rst Normal file
View File

@ -0,0 +1,9 @@
=================
Adding unit tests
=================
Unit tests in LBPM are implemented using ctest
* general overview
* launching unit tests for GPU (MPI flags etc.)

9
docs/source/examples/running.rst Normal file
View File

@ -0,0 +1,9 @@
============
Running LBPM
============
There are two main components to running LBPM simulators.
First is understanding how to launch MPI tasks on your system,
which depends on the particular implementation of MPI that you are using,
as well as other details of the local configuration. The second component is
understanding the LBPM input file structure.

25
docs/source/index.rst Normal file
View File

@ -0,0 +1,25 @@
.. LBPM documentation master file, created by
sphinx-quickstart on Thu May 20 12:19:14 2021.
You can adapt this file completely to your liking, but it should at least
contain the root `toctree` directive.
LBPM -- Documentation
===================================================
.. toctree::
:glob:
:maxdepth: 2
:caption: Contents:
install
examples/*
userGuide/*
developerGuide/*
publications/*
Indices and tables
==================
* :ref:`genindex`
* :ref:`modindex`
* :ref:`search`

10
docs/source/install.rst Normal file
View File

@ -0,0 +1,10 @@
============
Installation
============
LBPM requires CMake, MPI and HDF5 as required dependencies.
.. code-block:: bash
ls $LBPM_SOURCE/sample_scripts

13
docs/source/publications/publications.rst Normal file
View File

@ -0,0 +1,13 @@
============
Publications
============
* James E McClure, Zhe Li, Mark Berrill, Thomas Ramstad, "The LBPM software package for simulating multiphase flow on digital images of porous rocks" Computational Geosciences (25) 871–895 (2021) https://doi.org/10.1007/s10596-020-10028-9
* James E. McClure, Zhe Li, Adrian P. Sheppard, Cass T. Miller, "An adaptive volumetric flux boundary condition for lattice Boltzmann methods" Computers & Fluids (210) (2020) https://doi.org/10.1016/j.compfluid.2020.104670
* Y.D. Wang, T. Chung, R.T. Armstrong, J. McClure, T. Ramstad, P. Mostaghimi, "Accelerated Computation of Relative Permeability by Coupled Morphological and Direct Multiphase Flow Simulation" Journal of Computational Physics (401) (2020) https://doi.org/10.1016/j.jcp.2019.108966

13
docs/source/userGuide/IO/fileformat.rst Normal file
View File

@ -0,0 +1,13 @@
========================
I/O conventions for LBPM
========================
There are three main kinds of output file that are supported by LBPM.
* CSV files --
* formatted binary files --
* unformatted binary files --

5
docs/source/userGuide/analysis/analysis.rst Normal file
View File

@ -0,0 +1,5 @@
===========================
Internal Analysis Framework
===========================
placeholder for analysis

17
docs/source/userGuide/index.rst Normal file
View File

@ -0,0 +1,17 @@
###############################################################################
User Guide
###############################################################################
Welcome to the LBPM user guide.
.. toctree::
:glob:
:maxdepth: 2
models/*
analysis/*
visualization/*
IO/*

6
docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst Normal file
View File

@ -0,0 +1,6 @@
=============================================
Poisson-Boltzmann model
=============================================
The LBPM Poisson-Boltzmann solver is designed to solve the Poisson-Boltzmann equation
to solve for the electric field in an ionic fluid.

91
docs/source/userGuide/models/color/index.rst Normal file
View File

@ -0,0 +1,91 @@
###############################################################################
Color model
###############################################################################
The LBPM color model is implemented by combining a multi-relaxation time D3Q19
lattice Boltzmann equation (LBE) to solve for the momentum transport with two D3Q7
LBEs for the mass transport. The color model will obey strict mass and momentum
conservation while minimizing diffusive fluxes across the interface between fluids.
The color model is a good choice for modeling dense fluids that are strongly immiscible
(e.g. water-oil systems). Due to the strong anti-diffusion in the interface region,
the color model is not suitable for modeling processes such as Ostwald ripening that
depend on diffusive fluxes between fluid phases.
A typical command to launch the LBPM color simulator is as follows
```
mpirun -np $NUMPROCS lbpm_color_simulator input.db
```
where ``$NUMPROCS`` is the number of MPI processors to be used and ``input.db`` is
the name of the input database that provides the simulation parameters.
Note that the specific syntax to launch MPI tasks may vary depending on your system.
For additional details please refer to your local system documentation.
****************************
Simulation protocols
****************************
Simulation protocols are designed to make it simpler to design and execute common
computational experiments. Protocols will automatically determine boundary conditions
needed to perform a particular simulation. LBPM will internall set default simulation paramaters
that can be over-ridden to develop customized simulations.
.. toctree::
:glob:
:maxdepth: 2
protocols/*
***************************
Model parameters
***************************
The essential model parameters for the color model are
- :math:`\alpha` -- control the interfacial tension between fluids with key ``alpha``
- :math:`\beta` -- control the width of the interface with key ``beta``
- :math:`\tau_A` -- control the viscosity of fluid A with key ``tauA``
- :math:`\tau_B` -- control the viscosity of fluid B with key ``tauB``
****************************
Model Formulation
****************************
****************************
Boundary Conditions
****************************
The following external boundary conditions are supported by ``lbpm_color_simulator``
and can be set by setting the ``BC`` key values in the ``Domain`` section of the
input file database
- ``BC = 0`` -- fully periodic boundary conditions
- ``BC = 3`` -- constant pressure boundary condition
- ``BC = 4`` -- constant volumetric flux boundary condition
For ``BC = 0`` any mass that exits on one side of the domain will re-enter at the other
side. If the pore-structure for the image is tight, the mismatch between the inlet and
outlet can artificially reduce the permeability of the sample due to the blockage of
flow pathways at the boundary. LBPM includes an internal utility that will reduce the impact
of the boundary mismatch by eroding the solid labels within the inlet and outlet layers
(https://doi.org/10.1007/s10596-020-10028-9) to create a mixing layer.
The number mixing layers to use can be set using the key values in the ``Domain`` section
of the input database
- ``InletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the inlet
- ``OUtletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the outlet
For the other boundary conditions a thin reservoir of fluid (default ``3`` voxels)
is established at either side of the domain. The inlet is defined as the boundary face
where ``z = 0`` and the outlet is the boundary face where ``z = nprocz*nz``. By default a
reservoir of fluid A is established at the inlet and a reservoir of fluid B is established at
the outlet, each with a default thickness of three voxels. To over-ride the default label at
the inlet or outlet, the ``Domain`` section of the database may specify the following key values
- ``InletLayerPhase = 2`` -- establish a reservoir of component B at the inlet
- ``OutletLayerPhase = 1`` -- establish a reservoir of component A at the outlet

12
docs/source/userGuide/models/color/protocols/centrifuge.rst Normal file
View File

@ -0,0 +1,12 @@
======================================
Color model -- Centrifuge Protocol
======================================
The centrifuge protocol is designed to mimic SCAL centrifuge experiments that
are used to infer the capillary pressure. The centrifuge protocol
.. code-block:: bash
image_sequence = "image1.raw", "image2.raw"

63
docs/source/userGuide/models/color/protocols/coreFlooding.rst Normal file
View File

@ -0,0 +1,63 @@
======================================
Color model -- Core Flooding
======================================
The core flooding protocol is designed to mimic SCAL experiments where one
immiscible fluid is injected into the sample at a constant rate, displacing the
other fluid. The core flooding protocol relies on a flux boundary condition
to ensure that fluid is injected into the sample at a constant rate. The flux
boundary condition implements a time-varying pressure boundary condition that
adapts to ensure a constant volumetric flux. Details for the flux boundary
condition are available
(see: https://doi.org/10.1016/j.compfluid.2020.104670)
.. code-block:: bash
protocol = "core flooding"
To match experimental conditions, it is usually important to match the capillary
number, which is
.. math::
\mbox{Ca} = \frac{\mu u_z}{\gamma}
where :math:`\mu` is the dynamic viscosity, :math:`u_z` is the fluid
(usually water) velocity and :math:`\gamma` is the interfacial tension.
The volumetric flow rate is related to the fluid velocity based on
.. math::
Q_z = \epsilon C_{xy} u_z
where :math:`C_{xy}` is the cross-sectional area and :math:`\epsilon`
is the porosity. Given a particular experimental system
self-similar conditions can be determined for the lattice Boltzmann
simulation by matching the non-dimensional :math:`mbox{Ca}`. It is nearly
awlays advantageous for the timestep to be as large as possible so
that time-to-solution will be more favorable. This is accomplished by
* use a high value for the numerical surface tension (e.g. ``alpha=1.0e-2``)
* use a small value for the fluid viscosity (e.g. ``tau_w = 0.7`` and ``tau_n = 0.7`` )
* determine the volumetric flow rate needed to match :math:`\mbox{Ca}`
For the color LBM the interfacial tension is
:math:`\gamma = 6 \alpha` and the dynamic viscosity is :math:`\mu = \rho(\tau-1/2)/3`,
where the units are relative to the lattice spacing, timestep and mass
density. Agreemetn between the experimental and simulated values for
:math:`\mbox{Ca}` is ensured by setting the volumetric flux
.. math::
Q_z = \frac{\epsilon C_{xy} \gamma }{\mu} \mbox{Ca}
where the LB units of the volumetric flux will be voxels per timestep.
In some situations it may also be important to match other non-dimensional numbers,
such as the viscosity ratio, density ratio, and/or Ohnesorge/Laplace number. This
can be accomplished with an analogous procedure. Enforcing additional constraints
will necessarily restrict the LB parameters that can be used, which are ultimately
manifested as a constraint on the size of a timestep.

17
docs/source/userGuide/models/color/protocols/fractionalFlow.rst Normal file
View File

@ -0,0 +1,17 @@
==========================================
Color model -- Fractional Flow Protocol
==========================================
The fractional flow protocol is designed to perform steady-state relative
permeability simulations by using an internal routine to change the fluid
saturation by adding and subtracting mass to the fluid phases. The
mass density is updated for each fluid based on the locations where
the local rate of flow is high.
.. code-block:: bash
protocol = "fractional flow"

27
docs/source/userGuide/models/color/protocols/imageSequence.rst Normal file
View File

@ -0,0 +1,27 @@
======================================
Color model -- Image Sequence Protocol
======================================
The image sequence protocol is designed to perform a set steady-state
simulations based on a sequence of 3D (8-bit) images provided by the user.
The images might be the output of a previous LBPM simulation, a sequence of
(segmented) experimental data, or data generated from a custom routine.
The image sequence protocol will apply the same set of flow conditions
to all images in the sequence. This means
* the image labels and any associated properties are the same
* the external boundary conditions are the same
* the physical simulation parameters are the same
The image sequence protocol does not set boundary conditions by default.
It is up to the user to determine the flow condition, with the understanding
that the same set of will be applied to each image in the sequence.
To enable the image sequence protocol, the following keys should be set
within the ``Color`` section of the input database
.. code-block:: bash
protocol = "image sequence"
image_sequence = "image1.raw", "image2.raw"

18
docs/source/userGuide/models/color/protocols/shellAggregation.rst Normal file
View File

@ -0,0 +1,18 @@
==========================================
Color model -- Shell Aggregation Protocol
==========================================
The shell aggregation protocol is designed to perform steady-state relative
permeability simulations by using an internal routine to change the fluid
saturation by moving the interface. The basic design for the shell aggregation
protocol is described by Wang et al. ( https://doi.org/10.1016/j.jcp.2019.108966 ).
.. code-block:: bash
protocol = "shell aggregation"

6
docs/source/userGuide/models/freeEnergy/freeEnergy.rst Normal file
View File

@ -0,0 +1,6 @@
=============================================
Free energy model
=============================================
The LBPM free energy model is constructed to solve the Allen-Cahn equations,
which are typically used to describe liquid-gas systems.

7
docs/source/userGuide/models/greyscale/greyscale.rst Normal file
View File

@ -0,0 +1,7 @@
=============================================
Greyscale model model
=============================================
The LBPM greyscale lattice Boltzmann model is constructed to approximate the
solution of the Darcy-Brinkman equations in grey regions, coupled to a Navier-Stokes
solution in open regions.

23
docs/source/userGuide/models/index.rst Normal file
View File

@ -0,0 +1,23 @@
###############################################################################
LBPM model summary
###############################################################################
Currently supported lattice Boltzmann models
.. toctree::
:glob:
:maxdepth: 2
color/*
mrt/*
nernstPlanck/*
PoissonBoltzmann/*
greyscale/*
freeEnergy/*

6
docs/source/userGuide/models/mrt/mrt.rst Normal file
View File

@ -0,0 +1,6 @@
=============================================
MRT model
=============================================
The multi-relaxation time (MRT) lattice Boltzmann model is constructed to approximate the
solution of the Navier-Stokes equations.

6
docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst Normal file
View File

@ -0,0 +1,6 @@
=============================================
Nernst-Planck model
=============================================
The Nernst-Planck model is designed to model ion transport based on the
Nernst-Planck equation.

5
docs/source/userGuide/visualization/visit.rst Normal file
View File

@ -0,0 +1,5 @@
======================================
Visualizing simulation data with visit
======================================
placeholder for visit

928
models/ColorModel.cpp
View File

File diff suppressed because it is too large Load Diff

6
models/ColorModel.h
View File

@ -72,6 +72,8 @@ public:
double *ColorGrad;
double *Velocity;
double *Pressure;
void AssignComponentLabels(double *phase);
private:
Utilities::MPI comm;
@ -85,7 +87,6 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void AssignComponentLabels(double *phase);
double ImageInit(std::string filename);
double MorphInit(const double beta, const double morph_delta);
double SeedPhaseField(const double seed_water_in_oil);
@ -97,7 +98,10 @@ public:
FlowAdaptor(ScaLBL_ColorModel &M);
~FlowAdaptor();
double MoveInterface(ScaLBL_ColorModel &M);
double ImageInit(ScaLBL_ColorModel &M, std::string Filename);
double ShellAggregation(ScaLBL_ColorModel &M, const double delta_volume);
double UpdateFractionalFlow(ScaLBL_ColorModel &M);
double SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil);
void Flatten(ScaLBL_ColorModel &M);
DoubleArray phi;
DoubleArray phi_t;

73
models/FreeLeeModel.cpp
View File

@ -63,7 +63,6 @@ void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray
// Gets data (in optimized layout) from the HOST and stores in regular layout
// Primarly for debugging
int i,j,k,idx;
int n;
// initialize the array
regdata.fill(0.f);
@ -72,7 +71,6 @@ void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
n=k*Nx*Ny+j*Nx+i;
idx=Map(i,j,k);
if (!(idx<0)){
value=data[idx];
@ -414,8 +412,10 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels_ChemPotential_ColorGrad()
ERROR("Error: ComponentLabels and ComponentAffinity must be the same length! \n");
}
double label_count[NLABELS];
double label_count_global[NLABELS];
double *label_count;
double *label_count_global;
label_count = new double [NLABELS];
label_count_global = new double [NLABELS];
// Assign the labels
for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
@ -738,75 +738,10 @@ void ScaLBL_FreeLeeModel::Initialize_SingleFluid(){
if (Restart == true){
//TODO need to revise this function
//remove the phase-related part
// if (rank==0){
// printf("Reading restart file! \n");
// }
//
// // Read in the restart file to CPU buffers
// int *TmpMap;
// TmpMap = new int[Np];
//
// double *cPhi, *cDist, *cDen;
// cPhi = new double[N];
// cDen = new double[2*Np];
// cDist = new double[19*Np];
// ScaLBL_CopyToHost(TmpMap, dvcMap, Np*sizeof(int));
// //ScaLBL_CopyToHost(cPhi, Phi, N*sizeof(double));
//
// ifstream File(LocalRestartFile,ios::binary);
// int idx;
// double value,va,vb;
// for (int n=0; n<Np; n++){
// File.read((char*) &va, sizeof(va));
// File.read((char*) &vb, sizeof(vb));
// cDen[n] = va;
// cDen[Np+n] = vb;
// }
// for (int n=0; n<Np; n++){
// // Read the distributions
// for (int q=0; q<19; q++){
// File.read((char*) &value, sizeof(value));
// cDist[q*Np+n] = value;
// }
// }
// File.close();
//
// for (int n=0; n<ScaLBL_Comm->LastExterior(); n++){
// va = cDen[n];
// vb = cDen[Np + n];
// value = (va-vb)/(va+vb);
// idx = TmpMap[n];
// if (!(idx < 0) && idx<N)
// cPhi[idx] = value;
// }
// for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
// va = cDen[n];
// vb = cDen[Np + n];
// value = (va-vb)/(va+vb);
// idx = TmpMap[n];
// if (!(idx < 0) && idx<N)
// cPhi[idx] = value;
// }
//
// // Copy the restart data to the GPU
// ScaLBL_CopyToDevice(Den,cDen,2*Np*sizeof(double));
// ScaLBL_CopyToDevice(gqbar,cDist,19*Np*sizeof(double));
// ScaLBL_CopyToDevice(Phi,cPhi,N*sizeof(double));
// ScaLBL_Comm->Barrier();
// comm.barrier();
//
// if (rank==0) printf ("Initializing phase and density fields on device from Restart\n");
// //TODO the following function is to be updated.
// //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, 0, ScaLBL_Comm->LastExterior(), Np);
// //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
}
double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
int nprocs=nprocx*nprocy*nprocz;
int START_TIME = timestep;
int EXIT_TIME = min(returntime, timestepMax);

49
models/GreyscaleModel.cpp
View File

@ -150,7 +150,6 @@ void ScaLBL_GreyscaleModel::ReadInput(){
// Generate the signed distance map
// Initialize the domain and communication
Array<char> id_solid(Nx,Ny,Nz);
int count = 0;
// Solve for the position of the solid phase
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
@ -167,7 +166,6 @@ void ScaLBL_GreyscaleModel::ReadInput(){
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;
// Initialize distance to +/- 1
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
}
@ -199,11 +197,13 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
ERROR("Error: ComponentLabels and PorosityList must be the same length! \n");
}
double label_count[NLABELS];
double label_count_global[NLABELS];
// Assign the labels
for (int idx=0; idx<NLABELS; idx++) label_count[idx]=0;
double *label_count;
double *label_count_global;
label_count = new double [NLABELS];
label_count_global = new double [NLABELS];
for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
@ -211,7 +211,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
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++){
for (size_t idx=0; idx < NLABELS; idx++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
POROSITY=PorosityList[idx];
@ -242,7 +242,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
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++){
for (size_t idx=0; idx < NLABELS; idx++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
PERMEABILITY=PermeabilityList[idx];
@ -267,7 +267,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
// Set Dm to match Mask
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
for (int idx=0; idx<NLABELS; idx++) label_count_global[idx]=Dm->Comm.sumReduce( label_count[idx]);
for (size_t idx=0; idx<NLABELS; idx++) label_count_global[idx]=Dm->Comm.sumReduce( label_count[idx]);
//Initialize a weighted porosity after considering grey voxels
GreyPorosity=0.0;
for (unsigned int idx=0; idx<NLABELS; idx++){
@ -598,7 +598,7 @@ void ScaLBL_GreyscaleModel::Run(){
//double mass_loc,mass_glb;
//parameters for domain average
int64_t i,j,k,n,imin,jmin,kmin,kmax;
int64_t imin,jmin,kmin,kmax;
// If external boundary conditions are set, do not average over the inlet and outlet
kmin=1; kmax=Nz-1;
//In case user forgets to specify the inlet/outlet buffer layers for BC>0
@ -691,23 +691,22 @@ void ScaLBL_GreyscaleModel::Run(){
//double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag;
double absperm = h*h*mu*GreyPorosity*flow_rate / force_mag;
if (rank==0){
if (rank==0){
printf(" AbsPerm = %.5g [micron^2]\n",absperm);
bool WriteHeader=false;
FILE * log_file = fopen("Permeability.csv","r");
if (log_file != NULL)
fclose(log_file);
else
WriteHeader=true;
log_file = fopen("Permeability.csv","a");
if (WriteHeader)
fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n",
timestep,Fx,Fy,Fz,mu,h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz,absperm);
bool WriteHeader=false;
FILE * log_file = fopen("Permeability.csv","r");
if (log_file != NULL)
fclose(log_file);
else
WriteHeader=true;
log_file = fopen("Permeability.csv","a");
if (WriteHeader)
fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n");
fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu,
h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm);
fclose(log_file);
}
fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu,
h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm);
fclose(log_file);
}
}
if (timestep%visualization_interval==0){

78
models/IonModel.cpp
View File

@ -97,7 +97,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
}
else{
time_conv.clear();
for (int i=0; i<tau.size();i++){
for (size_t i=0; i<tau.size();i++){
time_conv.push_back((tau[i]-0.5)/k2_inv*(h*h*1.0e-12)/Di[i]);
}
}
@ -112,13 +112,13 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n");
}
else{
for (int i=0; i<IonDiffusivity.size();i++){
for (size_t i=0; i<IonDiffusivity.size();i++){
IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
}
}
}
else {
for (int i=0; i<IonDiffusivity.size();i++){
for (size_t i=0; i<IonDiffusivity.size();i++){
//convert ion diffusivity in physical unit to LB unit
IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
}
@ -141,13 +141,13 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
ERROR("Error: number_ion_species and IonConcentrationList must be the same length! \n");
}
else{
for (int i=0; i<IonConcentration.size();i++){
for (size_t i=0; i<IonConcentration.size();i++){
IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
else {
for (int i=0; i<IonConcentration.size();i++){
for (size_t i=0; i<IonConcentration.size();i++){
IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
@ -186,7 +186,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
else {
ERROR("Error: Non-periodic BCs are specified but InletValueList cannot be found! \n");
}
for (unsigned int i=0;i<BoundaryConditionInlet.size();i++){
for (size_t i=0;i<BoundaryConditionInlet.size();i++){
switch (BoundaryConditionInlet[i]){
case 1://fixed boundary ion concentration [mol/m^3]
Cin[i] = Cin[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
@ -220,7 +220,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
else {
ERROR("Error: Non-periodic BCs are specified but OutletValueList cannot be found! \n");
}
for (unsigned int i=0;i<BoundaryConditionOutlet.size();i++){
for (size_t i=0;i<BoundaryConditionOutlet.size();i++){
switch (BoundaryConditionOutlet[i]){
case 1://fixed boundary ion concentration [mol/m^3]
Cout[i] = Cout[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
@ -307,7 +307,7 @@ void ScaLBL_IonModel::ReadParams(string filename){
}
else{
time_conv.clear();
for (int i=0; i<tau.size();i++){
for (size_t i=0; i<tau.size();i++){
time_conv.push_back((tau[i]-0.5)/k2_inv*(h*h*1.0e-12)/Di[i]);
}
}
@ -322,13 +322,13 @@ void ScaLBL_IonModel::ReadParams(string filename){
ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n");
}
else{
for (int i=0; i<IonDiffusivity.size();i++){
for (size_t i=0; i<IonDiffusivity.size();i++){
IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
}
}
}
else {
for (int i=0; i<IonDiffusivity.size();i++){
for (size_t i=0; i<IonDiffusivity.size();i++){
//convert ion diffusivity in physical unit to LB unit
IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
}
@ -351,13 +351,13 @@ void ScaLBL_IonModel::ReadParams(string filename){
ERROR("Error: number_ion_species and IonConcentrationList must be the same length! \n");
}
else{
for (int i=0; i<IonConcentration.size();i++){
for (size_t i=0; i<IonConcentration.size();i++){
IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
else {
for (int i=0; i<IonConcentration.size();i++){
for (size_t i=0; i<IonConcentration.size();i++){
IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
@ -396,7 +396,7 @@ void ScaLBL_IonModel::ReadParams(string filename){
else {
ERROR("Error: Non-periodic BCs are specified but InletValueList cannot be found! \n");
}
for (unsigned int i=0;i<BoundaryConditionInlet.size();i++){
for (size_t i=0;i<BoundaryConditionInlet.size();i++){
switch (BoundaryConditionInlet[i]){
case 1://fixed boundary ion concentration [mol/m^3]
Cin[i] = Cin[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
@ -430,7 +430,7 @@ void ScaLBL_IonModel::ReadParams(string filename){
else {
ERROR("Error: Non-periodic BCs are specified but OutletValueList cannot be found! \n");
}
for (unsigned int i=0;i<BoundaryConditionOutlet.size();i++){
for (size_t i=0;i<BoundaryConditionOutlet.size();i++){
switch (BoundaryConditionOutlet[i]){
case 1://fixed boundary ion concentration [mol/m^3]
Cout[i] = Cout[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
@ -558,10 +558,12 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
ERROR("Error: LB Ion Solver: SolidLabels and SolidValues must be the same length! \n");
}
double label_count[NLABELS];
double label_count_global[NLABELS];
// Assign the labels
// Assign the labels
double *label_count;
double *label_count_global;
label_count = new double [NLABELS];
label_count_global = new double [NLABELS];
for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
for (int k=0;k<Nz;k++){
@ -571,7 +573,7 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
VALUE=Mask->id[n];
AFFINITY=0.f;
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
for (size_t idx=0; idx < NLABELS; idx++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
@ -701,23 +703,23 @@ void ScaLBL_IonModel::Initialize(){
auto File_ion = ion_db->getVector<std::string>( "IonConcentrationFile" );
double *Ci_host;
Ci_host = new double[number_ion_species*Np];
for (int ic=0; ic<number_ion_species; ic++){
for (size_t ic=0; ic<number_ion_species; ic++){
AssignIonConcentration_FromFile(&Ci_host[ic*Np],File_ion);
}
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species*sizeof(double)*Np);
comm.barrier();
for (int ic=0; ic<number_ion_species; ic++){
for (size_t ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic*Np*7],&Ci[ic*Np],Np);
}
delete [] Ci_host;
}
else{
for (int ic=0; ic<number_ion_species; ic++){
for (size_t ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init(&fq[ic*Np*7],&Ci[ic*Np],IonConcentration[ic],Np);
}
}
if (rank==0) printf ("LB Ion Solver: initializing charge density\n");
for (int ic=0; ic<number_ion_species; ic++){
for (size_t ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np);
}
@ -734,35 +736,35 @@ void ScaLBL_IonModel::Initialize(){
break;
}
for (int i=0; i<number_ion_species;i++){
for (size_t i=0; i<number_ion_species;i++){
switch (BoundaryConditionInlet[i]){
case 0:
if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %i is periodic \n",i+1);
if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %zu is periodic \n",i+1);
break;
case 1:
if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %i is concentration = %.5g [mol/m^3] \n",i+1,Cin[i]/(h*h*h*1.0e-18));
if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %zu is concentration = %.5g [mol/m^3] \n",i+1,Cin[i]/(h*h*h*1.0e-18));
break;
case 2:
if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %i is (inward) flux = %.5g [mol/m^2/sec] \n",i+1,Cin[i]/(h*h*1.0e-12)/time_conv[i]);
if (rank==0) printf("LB Ion Solver: inlet boundary for Ion %zu is (inward) flux = %.5g [mol/m^2/sec] \n",i+1,Cin[i]/(h*h*1.0e-12)/time_conv[i]);
break;
}
switch (BoundaryConditionOutlet[i]){
case 0:
if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %i is periodic \n",i+1);
if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %zu is periodic \n",i+1);
break;
case 1:
if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %i is concentration = %.5g [mol/m^3] \n",i+1,Cout[i]/(h*h*h*1.0e-18));
if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %zu is concentration = %.5g [mol/m^3] \n",i+1,Cout[i]/(h*h*h*1.0e-18));
break;
case 2:
if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %i is (inward) flux = %.5g [mol/m^2/sec] \n",i+1,Cout[i]/(h*h*1.0e-12)/time_conv[i]);
if (rank==0) printf("LB Ion Solver: outlet boundary for Ion %zu is (inward) flux = %.5g [mol/m^2/sec] \n",i+1,Cout[i]/(h*h*1.0e-12)/time_conv[i]);
break;
}
}
if (rank==0) printf("*****************************************************\n");
if (rank==0) printf("LB Ion Transport Solver: \n");
for (int i=0; i<number_ion_species;i++){
if (rank==0) printf(" Ion %i: LB relaxation tau = %.5g\n", i+1,tau[i]);
for (size_t i=0; i<number_ion_species;i++){
if (rank==0) printf(" Ion %zu: LB relaxation tau = %.5g\n", i+1,tau[i]);
if (rank==0) printf(" Time conversion factor: %.5g [sec/lt]\n", time_conv[i]);
if (rank==0) printf(" Internal iteration: %i [lt]\n", timestepMax[i]);
}
@ -777,7 +779,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//LB-related parameter
vector<double> rlx;
for (unsigned int ic=0;ic<tau.size();ic++){
for (size_t ic=0;ic<tau.size();ic++){
rlx.push_back(1.0/tau[ic]);
}
@ -786,7 +788,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//ScaLBL_Comm->Barrier(); comm.barrier();
//auto t1 = std::chrono::system_clock::now();
for (int ic=0; ic<number_ion_species; ic++){
for (size_t ic=0; ic<number_ion_species; ic++){
timestep=0;
while (timestep < timestepMax[ic]) {
//************************************************************************/
@ -881,7 +883,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
}
//Compute charge density for Poisson equation
for (int ic=0; ic<number_ion_species; ic++){
for (size_t ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np);
}
@ -901,7 +903,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//if (rank==0) printf("********************************************************\n");
}
void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const int ic){
void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const size_t ic){
//This function wirte out the data in a normal layout (by aggregating all decomposed domains)
ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration);
@ -913,13 +915,13 @@ void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const i
void ScaLBL_IonModel::getIonConcentration_debug(int timestep){
//This function write out decomposed data
DoubleArray PhaseField(Nx,Ny,Nz);
for (int ic=0; ic<number_ion_species; ic++){
for (size_t ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonConcentration_LB_to_Phys(PhaseField);
FILE *OUTFILE;
sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank);
sprintf(LocalRankFilename,"Ion%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE);
@ -986,7 +988,7 @@ double ScaLBL_IonModel::CalIonDenConvergence(vector<double> &ci_avg_previous){
Ci_host = new double[Np];
vector<double> error(number_ion_species,0.0);
for (int ic=0; ic<number_ion_species; ic++){
for (size_t ic=0; ic<number_ion_species; ic++){
ScaLBL_CopyToHost(Ci_host,&Ci[ic*Np],Np*sizeof(double));
double count_loc=0;

4
models/IonModel.h
View File

@ -34,7 +34,7 @@ public:
void Create();
void Initialize();
void Run(double *Velocity, double *ElectricField);
void getIonConcentration(DoubleArray &IonConcentration, const int ic);
void getIonConcentration(DoubleArray &IonConcentration, const size_t ic);
void getIonConcentration_debug(int timestep);
void DummyFluidVelocity();
void DummyElectricField();
@ -51,7 +51,7 @@ public:
double fluidVelx_dummy,fluidVely_dummy,fluidVelz_dummy;
double Ex_dummy,Ey_dummy,Ez_dummy;
int number_ion_species;
size_t number_ion_species;
vector<int> BoundaryConditionInlet;
vector<int> BoundaryConditionOutlet;
vector<double> IonDiffusivity;//User input unit [m^2/sec]

67
tests/DataAggregator.cpp
View File

@ -9,8 +9,8 @@
using namespace std;
int main(int argc, char **argv){
printf("Aggregating block data into single file \n");
printf("Aggregating block data into single file \n");
unsigned int Bx,By,Bz;
uint64_t Nx,Ny,Nz;
uint64_t N;
@ -22,25 +22,25 @@ int main(int argc, char **argv){
//Read the block size
if (argc>9){
printf("Input arguments accepted \n");
Nx = atoi(argv[1]);
Ny = atoi(argv[2]);
Nz = atoi(argv[3]);
x0 = atoi(argv[4]);
y0 = atoi(argv[5]);
z0 = atoi(argv[6]);
NX = atol(argv[7]);
NY = atol(argv[8]);
NZ = atol(argv[9]);
printf("Size %i X %i X %i \n",NX,NY,NZ);
fflush(stdout);
printf("Input arguments accepted \n");
Nx = atoi(argv[1]);
Ny = atoi(argv[2]);
Nz = atoi(argv[3]);
x0 = atoi(argv[4]);
y0 = atoi(argv[5]);
z0 = atoi(argv[6]);
NX = atol(argv[7]);
NY = atol(argv[8]);
NZ = atol(argv[9]);
printf("Size %llu X %llu X %llu \n",(unsigned long long) NX, (unsigned long long) NY, (unsigned long long) NZ);
fflush(stdout);
}
else{
printf("setting defaults \n");
Nx=Ny=Nz=1024;
x0=y0=z0=0;
NX=NY=8640;
NZ=6480;
printf("setting defaults \n");
Nx=Ny=Nz=1024;
x0=y0=z0=0;
NX=NY=8640;
NZ=6480;
}
//Bx = By = Bz = 9;
//Nx = Ny = Nz = 1024;
@ -53,29 +53,28 @@ int main(int argc, char **argv){
if (By>8) By=8;
if (Bz>8) Bz=8;
printf("System size (output) is: %i x %i x %i \n",NX,NY,NZ);
printf("Block size (read) is: %i x %i x %i \n",Nx,Ny,Nz);
printf("Starting location (read) is: %i, %i, %i \n", x0,y0,z0);
printf("Block number (read): %i x %i x %i \n",Bx,By,Bz);
printf("System size (output) is: %llu x %llu x %llu \n",(unsigned long long) NX,(unsigned long long) NY, (unsigned long long) NZ);
printf("Block size (read) is: %llu x %llu x %llu \n",(unsigned long long) Nx,(unsigned long long) Ny,(unsigned long long) Nz);
printf("Starting location (read) is: %llu, %llu, %llu \n", (unsigned long long) x0,(unsigned long long) y0,(unsigned long long) z0);
printf("Block number (read): %llu x %llu x %llu \n",(unsigned long long) Bx,(unsigned long long) By,(unsigned long long) Bz);
fflush(stdout);
// Filenames used
//char LocalRankString[8];
char LocalRankFilename[40];
char sx[2];
char sy[2];
char sz[2];
char tmpstr[10];
//sprintf(LocalRankString,"%05d",rank);
N = Nx*Ny*Nz;
N_full=NX*NY*NZ;
char *id;
id = new char [N];
char *ID;
ID = new char [N_full];
for (unsigned int bz=0; bz<Bz; bz++){
for (unsigned int by=0; by<By; by++){
for (unsigned int bx=0; bx<Bx; bx++){
@ -90,7 +89,7 @@ int main(int argc, char **argv){
FILE *IDFILE = fopen(LocalRankFilename,"rb");
readID=fread(id,1,N,IDFILE);
fclose(IDFILE);
printf("Loading data ... \n");
printf("Loading data: %zu bytes ... \n", readID);
// Unpack the data into the main array
for ( k=0;k<Nz;k++){
for ( j=0;j<Ny;j++){
@ -99,7 +98,7 @@ int main(int argc, char **argv){
y = by*Ny + j;
z = bz*Nz + k;
if ( x<NX && y<NY && z<NZ){
ID[z*NX*NY+y*NX+x] = id[k*Nx*Ny+j*Nx+i];
ID[z*NX*NY+y*NX+x] = id[k*Nx*Ny+j*Nx+i];
}
}
}
@ -111,11 +110,11 @@ int main(int argc, char **argv){
// Compute porosity
uint64_t count=0;
for (k=0; k<NZ; k++){
for (j=0; j<NY; j++){
for (i=0; i<NX; i++){
if (ID[k*NX*NY+j*NX+i] < 215) count++;
}
}
for (j=0; j<NY; j++){
for (i=0; i<NX; i++){
if (ID[k*NX*NY+j*NX+i] < 1) count++;
}
}
}
printf("Porosity is %f \n",double(count)/double(NX*NY*NZ));

14
tests/GenerateSphereTest.cpp
View File

@ -137,20 +137,18 @@ inline void MorphOpen(DoubleArray SignDist, char *id, Domain &Dm, int nx, int ny
int Nx = nx;
int Ny = ny;
int Nz = nz;
int imin,jmin,kmin,imax,jmax,kmax;
double sw_old=1.0;
double sw_new=1.0;
double sw_diff_old = 1.0;
double sw_diff_new = 1.0;
double sw_diff_old = 1.0;
double sw_diff_new = 1.0;
// Increase the critical radius until the target saturation is met
double deltaR=0.05; // amount to change the radius in voxel units
double Rcrit_old;
double Rcrit_new;
int imin,jmin,kmin,imax,jmax,kmax;
Rcrit_new = maxdistGlobal;
double Rcrit_new = maxdistGlobal;
double Rcrit_old = maxdistGlobal;
while (sw_new > SW)
{
sw_diff_old = sw_diff_new;

244
tests/lbpm_color_simulator.cpp
View File

@ -9,8 +9,6 @@
#include "common/Utilities.h"
#include "models/ColorModel.h"
//#define WRE_SURFACES
/*
* Simulator for two-phase flow in porous media
* James E. McClure 2013-2014
@ -24,100 +22,170 @@
int main( int argc, char **argv )
{
// Initialize
Utilities::startup( argc, argv );
// Initialize
Utilities::startup( argc, argv );
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
{ // 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";
}
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 );
NULL_USE( device );
ScaLBL_DeviceBarrier();
comm.barrier();
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 );
NULL_USE( device );
ScaLBL_DeviceBarrier();
comm.barrier();
PROFILE_ENABLE( 1 );
// PROFILE_ENABLE_TRACE();
// PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START( "Main" );
Utilities::setErrorHandlers();
PROFILE_ENABLE( 1 );
// PROFILE_ENABLE_TRACE();
// PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START( "Main" );
Utilities::setErrorHandlers();
auto filename = argv[1];
ScaLBL_ColorModel ColorModel( rank, nprocs, comm );
ColorModel.ReadParams( filename );
ColorModel.SetDomain();
ColorModel.ReadInput();
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
if (SimulationMode == "development"){
double MLUPS=0.0;
int timestep = 0;
/* flow adaptor keys to control */
auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" );
int MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
int SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 100000 );
auto filename = argv[1];
ScaLBL_ColorModel ColorModel( rank, nprocs, comm );
ColorModel.ReadParams( filename );
ColorModel.SetDomain();
ColorModel.ReadInput();
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
int ANALYSIS_INTERVAL = ColorModel.timestepMax;
if (ColorModel.analysis_db->keyExists( "analysis_interval" )){
ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
}
/* Launch the simulation */
FlowAdaptor Adapt(ColorModel);
runAnalysis analysis(ColorModel);
while (ColorModel.timestep < ColorModel.timestepMax){
/* this will run steady points */
timestep += MAX_STEADY_TIME;
MLUPS = ColorModel.Run(timestep);
if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
Adapt.UpdateFractionalFlow(ColorModel);
/* apply timestep skipping algorithm to accelerate steady-state */
int skip_time = 0;
timestep = ColorModel.timestep;
while (skip_time < SKIP_TIMESTEPS){
timestep += ANALYSIS_INTERVAL;
MLUPS = ColorModel.Run(timestep);
Adapt.MoveInterface(ColorModel);
skip_time += ANALYSIS_INTERVAL;
}
//Adapt.Flatten(ColorModel);
if (SimulationMode == "development"){
double MLUPS=0.0;
int timestep = 0;
bool ContinueSimulation = true;
/* Variables for simulation protocols */
auto PROTOCOL = ColorModel.color_db->getWithDefault<std::string>( "protocol", "none" );
/* image sequence protocol */
int IMAGE_INDEX = 0;
int IMAGE_COUNT = 0;
std::vector<std::string> ImageList;
/* flow adaptor keys to control behavior */
int SKIP_TIMESTEPS = 0;
int MAX_STEADY_TIME = 1000000;
double ENDPOINT_THRESHOLD = 0.1;
double FRACTIONAL_FLOW_INCREMENT = 0.0; // this will skip the flow adaptor if not enabled
double SEED_WATER = 0.0;
if (ColorModel.db->keyExists( "FlowAdaptor" )){
auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" );
MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 50000 );
ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
/* protocol specific key values */
if (PROTOCOL == "fractional flow")
FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
if (PROTOCOL == "seed water")
SEED_WATER = flow_db->getWithDefault<double>( "seed_water", 0.01);
}
/* analysis keys*/
int ANALYSIS_INTERVAL = ColorModel.timestepMax;
if (ColorModel.analysis_db->keyExists( "analysis_interval" )){
ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
}
/* Launch the simulation */
FlowAdaptor Adapt(ColorModel);
runAnalysis analysis(ColorModel);
while (ContinueSimulation){
/* this will run steady points */
timestep += MAX_STEADY_TIME;
MLUPS = ColorModel.Run(timestep);
if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
if (ColorModel.timestep > ColorModel.timestepMax){
ContinueSimulation = false;
}
/* Load a new image if image sequence is specified */
if (PROTOCOL == "image sequence"){
IMAGE_INDEX++;
if (IMAGE_INDEX < IMAGE_COUNT){
std::string next_image = ImageList[IMAGE_INDEX];
if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX);
ColorModel.color_db->putScalar<int>("image_index",IMAGE_INDEX);
Adapt.ImageInit(ColorModel, next_image);
}
else{
if (rank==0) printf("Finished simulating image sequence \n");
ColorModel.timestep = ColorModel.timestepMax;
}
}
/*********************************************************/
/* update the fluid configuration with the flow adapter */
int skip_time = 0;
timestep = ColorModel.timestep;
/* get the averaged flow measures computed internally for the last simulation point*/
double SaturationChange = 0.0;
double volB = ColorModel.Averages->gwb.V;
double volA = ColorModel.Averages->gnb.V;
double initialSaturation = volB/(volA + volB);
double vA_x = ColorModel.Averages->gnb.Px/ColorModel.Averages->gnb.M;
double vA_y = ColorModel.Averages->gnb.Py/ColorModel.Averages->gnb.M;
double vA_z = ColorModel.Averages->gnb.Pz/ColorModel.Averages->gnb.M;
double vB_x = ColorModel.Averages->gwb.Px/ColorModel.Averages->gwb.M;
double vB_y = ColorModel.Averages->gwb.Py/ColorModel.Averages->gwb.M;
double vB_z = ColorModel.Averages->gwb.Pz/ColorModel.Averages->gwb.M;
double speedA = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
/* stop simulation if previous point was sufficiently close to the endpoint*/
if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false;
if (ContinueSimulation){
while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
timestep += ANALYSIS_INTERVAL;
if (PROTOCOL == "fractional flow") {
Adapt.UpdateFractionalFlow(ColorModel);
}
else if (PROTOCOL == "shell aggregation"){
double target_volume_change = FRACTIONAL_FLOW_INCREMENT*initialSaturation - SaturationChange;
Adapt.ShellAggregation(ColorModel,target_volume_change);
}
else if (PROTOCOL == "seed water"){
Adapt.SeedPhaseField(ColorModel,SEED_WATER);
}
/* Run some LBM timesteps to let the system relax a bit */
MLUPS = ColorModel.Run(timestep);
/* Recompute the volume fraction now that the system has adjusted */
double volB = ColorModel.Averages->gwb.V;
double volA = ColorModel.Averages->gnb.V;
SaturationChange = volB/(volA + volB) - initialSaturation;
skip_time += ANALYSIS_INTERVAL;
}
if (rank==0) printf(" ********************************************************************* \n");
if (rank==0) printf(" Updated fractional flow with saturation change = %f \n", SaturationChange);
if (rank==0) printf(" Used protocol = %s \n", PROTOCOL.c_str());
if (rank==0) printf(" ********************************************************************* \n");
}
/*********************************************************/
}
}
else
ColorModel.Run();
}
ColorModel.WriteDebug();
}
else
ColorModel.Run();
PROFILE_STOP( "Main" );
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
auto level = db->getWithDefault<int>( "TimerLevel", 1 );
PROFILE_SAVE( file, level );
// ****************************************************
PROFILE_STOP( "Main" );
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
auto level = db->getWithDefault<int>( "TimerLevel", 1 );
NULL_USE(level);
PROFILE_SAVE( file, level );
// ****************************************************
} // Limit scope so variables that contain communicators will free before MPI_Finialize
} // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::shutdown();
return 0;
Utilities::shutdown();
return 0;
}

1
tests/lbpm_freelee_SingleFluidBGK_simulator.cpp
View File

@ -59,6 +59,7 @@ int main( int argc, char **argv )
PROFILE_STOP("Main");
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_freelee_SingleFluidBGK_simulator" );
auto level = db->getWithDefault<int>( "TimerLevel", 1 );
NULL_USE(level);
PROFILE_SAVE( file,level );
// ****************************************************

1
tests/lbpm_freelee_simulator.cpp
View File

@ -83,6 +83,7 @@ int main( int argc, char **argv )
PROFILE_STOP("Main");
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_freelee_simulator" );
auto level = db->getWithDefault<int>( "TimerLevel", 1 );
NULL_USE(level);
PROFILE_SAVE( file,level );
// ****************************************************