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

Merge branch 'master' into morph_refactor

Browse Source
This commit is contained in:
James McClure 2021-10-27 09:03:13 -04:00
commit 7ed5be99ba
52 changed files with 2076 additions and 2727 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

255
analysis/ElectroChemistry.cpp
View File

@ -10,13 +10,25 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr <Domain> dm):
ChemicalPotential.resize(Nx,Ny,Nz); ChemicalPotential.fill(0);
ElectricalPotential.resize(Nx,Ny,Nz); ElectricalPotential.fill(0);
ElectricalField_x.resize(Nx,Ny,Nz); ElectricalField_x.fill(0);
ElectricalField_y.resize(Nx,Ny,Nz); ElectricalField_y.fill(0);
ElectricalField_z.resize(Nx,Ny,Nz); ElectricalField_z.fill(0);
Pressure.resize(Nx,Ny,Nz); Pressure.fill(0);
Rho.resize(Nx,Ny,Nz); Rho.fill(0);
Vel_x.resize(Nx,Ny,Nz); Vel_x.fill(0); // Gradient of the phase indicator field
Vel_y.resize(Nx,Ny,Nz); Vel_y.fill(0);
Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0);
SDs.resize(Nx,Ny,Nz); SDs.fill(0);
IonFluxDiffusive_x.resize(Nx,Ny,Nz); IonFluxDiffusive_x.fill(0);
IonFluxDiffusive_y.resize(Nx,Ny,Nz); IonFluxDiffusive_y.fill(0);
IonFluxDiffusive_z.resize(Nx,Ny,Nz); IonFluxDiffusive_z.fill(0);
IonFluxAdvective_x.resize(Nx,Ny,Nz); IonFluxAdvective_x.fill(0);
IonFluxAdvective_y.resize(Nx,Ny,Nz); IonFluxAdvective_y.fill(0);
IonFluxAdvective_z.resize(Nx,Ny,Nz); IonFluxAdvective_z.fill(0);
IonFluxElectrical_x.resize(Nx,Ny,Nz); IonFluxElectrical_x.fill(0);
IonFluxElectrical_y.resize(Nx,Ny,Nz); IonFluxElectrical_y.fill(0);
IonFluxElectrical_z.resize(Nx,Ny,Nz); IonFluxElectrical_z.fill(0);
if (Dm->rank()==0){
bool WriteHeader=false;
TIMELOG = fopen("electrokinetic.csv","r");
@ -159,21 +171,55 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
visData[0].meshName = "domain";
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>();
//electric potential
auto ElectricPotentialVar = std::make_shared<IO::Variable>();
//electric field
auto ElectricFieldVar_x = std::make_shared<IO::Variable>();
auto ElectricFieldVar_y = std::make_shared<IO::Variable>();
auto ElectricFieldVar_z = std::make_shared<IO::Variable>();
//ion concentration
std::vector<shared_ptr<IO::Variable>> IonConcentration;
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
IonConcentration.push_back(std::make_shared<IO::Variable>());
}
//fluid velocity
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
// diffusive ion flux
std::vector<shared_ptr<IO::Variable>> IonFluxDiffusive;
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
//push in x-,y-, and z-component for each ion species
IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
IonFluxDiffusive.push_back(std::make_shared<IO::Variable>());
}
// advective ion flux
std::vector<shared_ptr<IO::Variable>> IonFluxAdvective;
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
//push in x-,y-, and z-component for each ion species
IonFluxAdvective.push_back(std::make_shared<IO::Variable>());
IonFluxAdvective.push_back(std::make_shared<IO::Variable>());
IonFluxAdvective.push_back(std::make_shared<IO::Variable>());
}
// electro-migrational ion flux
std::vector<shared_ptr<IO::Variable>> IonFluxElectrical;
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
//push in x-,y-, and z-component for each ion species
IonFluxElectrical.push_back(std::make_shared<IO::Variable>());
IonFluxElectrical.push_back(std::make_shared<IO::Variable>());
IonFluxElectrical.push_back(std::make_shared<IO::Variable>());
}
//--------------------------------------------------------------------------------------------------------------------
//-------------------------------------Create Names for Variables------------------------------------------------------
if (vis_db->getWithDefault<bool>( "save_electric_potential", true )){
ElectricPotential->name = "ElectricPotential";
ElectricPotential->type = IO::VariableType::VolumeVariable;
ElectricPotential->dim = 1;
ElectricPotential->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricPotential);
ElectricPotentialVar->name = "ElectricPotential";
ElectricPotentialVar->type = IO::VariableType::VolumeVariable;
ElectricPotentialVar->dim = 1;
ElectricPotentialVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricPotentialVar);
}
if (vis_db->getWithDefault<bool>( "save_concentration", true )){
@ -204,7 +250,105 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VzVar);
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_diffusive", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxDiffusive_x",ion+1);
IonFluxDiffusive[3*ion+0]->name = VisName;
IonFluxDiffusive[3*ion+0]->type = IO::VariableType::VolumeVariable;
IonFluxDiffusive[3*ion+0]->dim = 1;
IonFluxDiffusive[3*ion+0]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxDiffusive[3*ion+0]);
// y-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxDiffusive_y",ion+1);
IonFluxDiffusive[3*ion+1]->name = VisName;
IonFluxDiffusive[3*ion+1]->type = IO::VariableType::VolumeVariable;
IonFluxDiffusive[3*ion+1]->dim = 1;
IonFluxDiffusive[3*ion+1]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxDiffusive[3*ion+1]);
// z-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxDiffusive_z",ion+1);
IonFluxDiffusive[3*ion+2]->name = VisName;
IonFluxDiffusive[3*ion+2]->type = IO::VariableType::VolumeVariable;
IonFluxDiffusive[3*ion+2]->dim = 1;
IonFluxDiffusive[3*ion+2]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxDiffusive[3*ion+2]);
}
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_advective", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of advective flux
sprintf(VisName,"Ion%zu_FluxAdvective_x",ion+1);
IonFluxAdvective[3*ion+0]->name = VisName;
IonFluxAdvective[3*ion+0]->type = IO::VariableType::VolumeVariable;
IonFluxAdvective[3*ion+0]->dim = 1;
IonFluxAdvective[3*ion+0]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxAdvective[3*ion+0]);
// y-component of advective flux
sprintf(VisName,"Ion%zu_FluxAdvective_y",ion+1);
IonFluxAdvective[3*ion+1]->name = VisName;
IonFluxAdvective[3*ion+1]->type = IO::VariableType::VolumeVariable;
IonFluxAdvective[3*ion+1]->dim = 1;
IonFluxAdvective[3*ion+1]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxAdvective[3*ion+1]);
// z-component of advective flux
sprintf(VisName,"Ion%zu_FluxAdvective_z",ion+1);
IonFluxAdvective[3*ion+2]->name = VisName;
IonFluxAdvective[3*ion+2]->type = IO::VariableType::VolumeVariable;
IonFluxAdvective[3*ion+2]->dim = 1;
IonFluxAdvective[3*ion+2]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxAdvective[3*ion+2]);
}
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_electrical", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of electro-migrational flux
sprintf(VisName,"Ion%zu_FluxElectrical_x",ion+1);
IonFluxElectrical[3*ion+0]->name = VisName;
IonFluxElectrical[3*ion+0]->type = IO::VariableType::VolumeVariable;
IonFluxElectrical[3*ion+0]->dim = 1;
IonFluxElectrical[3*ion+0]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxElectrical[3*ion+0]);
// y-component of electro-migrational flux
sprintf(VisName,"Ion%zu_FluxElectrical_y",ion+1);
IonFluxElectrical[3*ion+1]->name = VisName;
IonFluxElectrical[3*ion+1]->type = IO::VariableType::VolumeVariable;
IonFluxElectrical[3*ion+1]->dim = 1;
IonFluxElectrical[3*ion+1]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxElectrical[3*ion+1]);
// z-component of electro-migrational flux
sprintf(VisName,"Ion%zu_FluxElectrical_z",ion+1);
IonFluxElectrical[3*ion+2]->name = VisName;
IonFluxElectrical[3*ion+2]->type = IO::VariableType::VolumeVariable;
IonFluxElectrical[3*ion+2]->dim = 1;
IonFluxElectrical[3*ion+2]->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(IonFluxElectrical[3*ion+2]);
}
}
if (vis_db->getWithDefault<bool>( "save_electric_field", false )){
ElectricFieldVar_x->name = "ElectricField_x";
ElectricFieldVar_x->type = IO::VariableType::VolumeVariable;
ElectricFieldVar_x->dim = 1;
ElectricFieldVar_x->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricFieldVar_x);
ElectricFieldVar_y->name = "ElectricField_y";
ElectricFieldVar_y->type = IO::VariableType::VolumeVariable;
ElectricFieldVar_y->dim = 1;
ElectricFieldVar_y->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricFieldVar_y);
ElectricFieldVar_z->name = "ElectricField_z";
ElectricFieldVar_z->type = IO::VariableType::VolumeVariable;
ElectricFieldVar_z->dim = 1;
ElectricFieldVar_z->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(ElectricFieldVar_z);
}
//--------------------------------------------------------------------------------------------------------------------
//------------------------------------Save All Variables--------------------------------------------------------------
if (vis_db->getWithDefault<bool>( "save_electric_potential", true )){
ASSERT(visData[0].vars[0]->name=="ElectricPotential");
Poisson.getElectricPotential(ElectricalPotential);
@ -215,7 +359,7 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
if (vis_db->getWithDefault<bool>( "save_concentration", true )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
sprintf(VisName,"IonConcentration_%zu",ion+1);
IonConcentration[ion]->name = VisName;
//IonConcentration[ion]->name = VisName;
ASSERT(visData[0].vars[1+ion]->name==VisName);
Array<double>& IonConcentrationData = visData[0].vars[1+ion]->data;
Ion.getIonConcentration(Rho,ion);
@ -235,10 +379,101 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P
fillData.copy(Vel_y,VelyData);
fillData.copy(Vel_z,VelzData);
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_diffusive", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxDiffusive_x",ion+1);
//IonFluxDiffusive[3*ion+0]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species+3*ion+0]->name==VisName);
// y-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxDiffusive_y",ion+1);
//IonFluxDiffusive[3*ion+1]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species+3*ion+1]->name==VisName);
// z-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxDiffusive_z",ion+1);
//IonFluxDiffusive[3*ion+2]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species+3*ion+2]->name==VisName);
Array<double>& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species+3*ion+0]->data;
Array<double>& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species+3*ion+1]->data;
Array<double>& IonFluxData_z = visData[0].vars[4+Ion.number_ion_species+3*ion+2]->data;
Ion.getIonFluxDiffusive(IonFluxDiffusive_x,IonFluxDiffusive_y,IonFluxDiffusive_z,ion);
fillData.copy(IonFluxDiffusive_x,IonFluxData_x);
fillData.copy(IonFluxDiffusive_y,IonFluxData_y);
fillData.copy(IonFluxDiffusive_z,IonFluxData_z);
}
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_advective", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxAdvective_x",ion+1);
//IonFluxDiffusive[3*ion+0]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+0]->name==VisName);
// y-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxAdvective_y",ion+1);
//IonFluxDiffusive[3*ion+1]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+1]->name==VisName);
// z-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxAdvective_z",ion+1);
//IonFluxDiffusive[3*ion+2]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+2]->name==VisName);
Array<double>& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+0]->data;
Array<double>& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+1]->data;
Array<double>& IonFluxData_z = visData[0].vars[4+Ion.number_ion_species*(1+3)+3*ion+2]->data;
Ion.getIonFluxAdvective(IonFluxAdvective_x,IonFluxAdvective_y,IonFluxAdvective_z,ion);
fillData.copy(IonFluxAdvective_x,IonFluxData_x);
fillData.copy(IonFluxAdvective_y,IonFluxData_y);
fillData.copy(IonFluxAdvective_z,IonFluxData_z);
}
}
if (vis_db->getWithDefault<bool>( "save_ion_flux_electrical", false )){
for (size_t ion=0; ion<Ion.number_ion_species; ion++){
// x-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxElectrical_x",ion+1);
//IonFluxDiffusive[3*ion+0]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+0]->name==VisName);
// y-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxElectrical_y",ion+1);
//IonFluxDiffusive[3*ion+1]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+1]->name==VisName);
// z-component of diffusive flux
sprintf(VisName,"Ion%zu_FluxElectrical_z",ion+1);
//IonFluxDiffusive[3*ion+2]->name = VisName;
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+2]->name==VisName);
Array<double>& IonFluxData_x = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+0]->data;
Array<double>& IonFluxData_y = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+1]->data;
Array<double>& IonFluxData_z = visData[0].vars[4+Ion.number_ion_species*(1+6)+3*ion+2]->data;
Ion.getIonFluxElectrical(IonFluxElectrical_x,IonFluxElectrical_y,IonFluxElectrical_z,ion);
fillData.copy(IonFluxElectrical_x,IonFluxData_x);
fillData.copy(IonFluxElectrical_y,IonFluxData_y);
fillData.copy(IonFluxElectrical_z,IonFluxData_z);
}
}
if (vis_db->getWithDefault<bool>( "save_electric_field", false )){
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+9)+0]->name=="ElectricField_x");
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+9)+1]->name=="ElectricField_y");
ASSERT(visData[0].vars[4+Ion.number_ion_species*(1+9)+2]->name=="ElectricField_z");
Poisson.getElectricField(ElectricalField_x, ElectricalField_y, ElectricalField_z);
Array<double>& ElectricalFieldxData = visData[0].vars[4+Ion.number_ion_species*(1+9)+0]->data;
Array<double>& ElectricalFieldyData = visData[0].vars[4+Ion.number_ion_species*(1+9)+1]->data;
Array<double>& ElectricalFieldzData = visData[0].vars[4+Ion.number_ion_species*(1+9)+2]->data;
fillData.copy(ElectricalField_x,ElectricalFieldxData);
fillData.copy(ElectricalField_y,ElectricalFieldyData);
fillData.copy(ElectricalField_z,ElectricalFieldzData);
}
if (vis_db->getWithDefault<bool>( "write_silo", true ))
IO::writeData( timestep, visData, Dm->Comm );
//--------------------------------------------------------------------------------------------------------------------
/* if (vis_db->getWithDefault<bool>( "save_8bit_raw", true )){
char CurrentIDFilename[40];
sprintf(CurrentIDFilename,"id_t%d.raw",timestep);

12
analysis/ElectroChemistry.h
View File

@ -35,11 +35,23 @@ public:
DoubleArray Rho; // density field
DoubleArray ChemicalPotential; // density field
DoubleArray ElectricalPotential; // density field
DoubleArray ElectricalField_x; // density field
DoubleArray ElectricalField_y; // density field
DoubleArray ElectricalField_z; // density field
DoubleArray Pressure; // pressure field
DoubleArray Vel_x; // velocity field
DoubleArray Vel_y;
DoubleArray Vel_z;
DoubleArray SDs;
DoubleArray IonFluxDiffusive_x; //ion diffusive flux components
DoubleArray IonFluxDiffusive_y;
DoubleArray IonFluxDiffusive_z;
DoubleArray IonFluxAdvective_x; //ion advective flux components
DoubleArray IonFluxAdvective_y;
DoubleArray IonFluxAdvective_z;
DoubleArray IonFluxElectrical_x; //ion electromigration flux components
DoubleArray IonFluxElectrical_y;
DoubleArray IonFluxElectrical_z;
ElectroChemistryAnalyzer(std::shared_ptr <Domain> Dm);
~ElectroChemistryAnalyzer();

6
analysis/FlowAdaptor.cpp
View File

@ -68,11 +68,11 @@ double FlowAdaptor::ImageInit(ScaLBL_ColorModel &M, std::string Filename){
double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
double MASS_FRACTION_CHANGE = 0.01;
double MASS_FRACTION_CHANGE = 0.006;
double FRACTIONAL_FLOW_EPSILON = 5e-6;
if (M.db->keyExists( "FlowAdaptor" )){
auto flow_db = M.db->getDatabase( "FlowAdaptor" );
MASS_FRACTION_CHANGE = flow_db->getWithDefault<double>( "mass_fraction_factor", 0.01);
MASS_FRACTION_CHANGE = flow_db->getWithDefault<double>( "mass_fraction_factor", 0.006);
FRACTIONAL_FLOW_EPSILON = flow_db->getWithDefault<double>( "fractional_flow_epsilon", 5e-6);
}
int Np = M.Np;
@ -511,4 +511,4 @@ double FlowAdaptor::SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water
ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
return(mass_loss);
}
}

20
analysis/GreyPhase.cpp
View File

@ -19,6 +19,7 @@ GreyPhaseAnalysis::GreyPhaseAnalysis(std::shared_ptr <Domain> dm):
Vel_x.resize(Nx,Ny,Nz); Vel_x.fill(0); // Gradient of the phase indicator field
Vel_y.resize(Nx,Ny,Nz); Vel_y.fill(0);
Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0);
MobilityRatio.resize(Nx,Ny,Nz); MobilityRatio.fill(0);
//.........................................
if (Dm->rank()==0){
@ -89,14 +90,17 @@ void GreyPhaseAnalysis::Basic(){
double nB = Rho_w(n);
double phi = (nA-nB)/(nA+nB);
double porosity = Porosity(n);
Water_local.M += rho_w*nB*porosity;
Water_local.Px += porosity*rho_w*nB*Vel_x(n);
Water_local.Py += porosity*rho_w*nB*Vel_y(n);
Water_local.Pz += porosity*rho_w*nB*Vel_z(n);
Oil_local.M += rho_n*nA*porosity;
Oil_local.Px += porosity*rho_n*nA*Vel_x(n);
Oil_local.Py += porosity*rho_n*nA*Vel_y(n);
Oil_local.Pz += porosity*rho_n*nA*Vel_z(n);
double mobility_ratio = MobilityRatio(n);
Water_local.M += nB*porosity;
Water_local.Px += porosity*(nA+nB)*Vel_x(n)*0.5*(1.0-mobility_ratio);
Water_local.Py += porosity*(nA+nB)*Vel_y(n)*0.5*(1.0-mobility_ratio);
Water_local.Pz += porosity*(nA+nB)*Vel_z(n)*0.5*(1.0-mobility_ratio);
Oil_local.M += nA*porosity;
Oil_local.Px += porosity*(nA+nB)*Vel_x(n)*0.5*(1.0+mobility_ratio);
Oil_local.Py += porosity*(nA+nB)*Vel_y(n)*0.5*(1.0+mobility_ratio);
Oil_local.Pz += porosity*(nA+nB)*Vel_z(n)*0.5*(1.0+mobility_ratio);
if ( phi > 0.99 ){
Oil_local.p += Pressure(n);

1
analysis/GreyPhase.h
View File

@ -71,6 +71,7 @@ public:
DoubleArray Vel_x; // velocity field
DoubleArray Vel_y;
DoubleArray Vel_z;
DoubleArray MobilityRatio;
GreyPhaseAnalysis(std::shared_ptr <Domain> Dm);
~GreyPhaseAnalysis();

200
common/Domain.h
View File

@ -16,87 +16,98 @@
#include "common/MPI.h"
#include "common/Communication.h"
#include "common/Database.h"
class Domain;
template<class TYPE> class PatchData;
/**
* @file Domain.h
* \brief Parallel Domain data structures and helper functions
*/
//! Class to hold information about a box
/**
* \class Box
*
* @details
* information about a box
*/
class Box {
public:
int ifirst[3];
int ilast[3];
};
class Patch;
enum class DataLocation { CPU, DEVICE };
/**
* \class Domain
*
* @details
* the Domain class includes basic information to distribution 3D image data to multiple processes using MPI.
* A regular domain decomposision is performed, with each MPI process getting a [Nx,Ny,Nz] sub-domain.
* 8-bit image labels are retained internally.
* The domain class resides on the CPU and provides utilities to support CPU-based analysis.
* GPU-based data structures should be constructed separately but may utilize information that the Domain class provides.
*/
//! Class to hold information about a patch
class Patch {
public:
//! Empty constructor
Patch() = delete;
//! Copy constructor
Patch( const Patch& ) = delete;
//! Assignment operator
Patch& operator=( const Patch& ) = delete;
//! Return the box for the patch
inline const Box& getBox() const { return d_box; }
//! Create patch data
template<class TYPE>
std::shared_ptr<PatchData<TYPE>> createPatchData( DataLocation location ) const;
private:
Box d_box;
int d_owner;
Domain *d_domain;
};
//! Class to hold domain info
class Domain{
public:
//! Default constructor
/**
* \brief Constructor
* @param db input database
* @param Communicator MPI communicator
*/
Domain( std::shared_ptr<Database> db, const Utilities::MPI& Communicator);
//! Obsolete constructor
/**
* \brief Obsolete constructor
*/
Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
double lx, double ly, double lz, int BC);
//! Empty constructor
/**
* \brief Empty constructor
*/
Domain() = delete;
//! Copy constructor
/**
* \brief Copy constructor
*/
Domain( const Domain& ) = delete;
//! Assignment operator
/**
* \brief Assignment operator
*/
Domain& operator=( const Domain& ) = delete;
//! Destructor
/**
* \brief Destructor
*/
~Domain();
//! Get the database
/**
* \brief Get the database
*/
inline std::shared_ptr<const Database> getDatabase() const { return d_db; }
//! Get the domain box
/**
* \brief Get the domain box
*/
inline const Box& getBox() const { return d_box; }
//! Get local patch
/**
* \brief Get local patch
*/
inline const Patch& getLocalPatch() const { return *d_localPatch; }
//! Get all patches
/**
* \brief Get all patches
*/
inline const std::vector<Patch>& getAllPatch() const { return d_patches; }
private:
/**
* \brief initialize from database
*/
void initialize( std::shared_ptr<Database> db );
std::shared_ptr<Database> d_db;
@ -124,6 +135,9 @@ public: // Public variables (need to create accessors instead)
//**********************************
// MPI ranks for all 18 neighbors
//**********************************
/**
* \brief Compute the porosity based on the current domain id file
*/
inline double Porosity() const { return porosity; }
inline int iproc() const { return rank_info.ix; }
inline int jproc() const { return rank_info.jy; }
@ -165,22 +179,78 @@ public: // Public variables (need to create accessors instead)
// Solid indicator function
std::vector<signed char> id;
/**
* \brief Read domain IDs from file
*/
void ReadIDs();
/**
* \brief Compute the porosity
*/
void ComputePorosity();
/**
* \brief Read image and perform domain decomposition
* @param filename - name of file to read IDs
*/
void Decomp( const std::string& filename );
/**
* \brief Perform a halo exchange using MPI
* @param Mesh - array data that holds scalar values
*/
void CommunicateMeshHalo(DoubleArray &Mesh);
/**
* \brief Initialize communication data structures within Domain object.
* This routine needs to be called before the communication functionality will work
*/
void CommInit();
/**
* \brief Count number of pore nodes (labels > 1)
*/
int PoreCount();
/**
* \brief Read array data from a file and distribute to local arrays for each MPI process
* @param Filename - name of the file to read the data
* @param Datatype - data type to use
* @param UserData - Array to store the values that are read
*/
void ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData);
/**
* \brief Aggregate labels from all MPI processes and write to a file
* @param filename - name of the file to write
*/
void AggregateLabels( const std::string& filename );
/**
* \brief Aggregate user provided array from all MPI processes and write to a single file
* @param filename - name of the file to write
* @param UserData - array data to aggregate and write
*/
void AggregateLabels( const std::string& filename, DoubleArray &UserData );
private:
/**
* \brief Pack halo data for 8-bit integer
* @param list - list of values in the halo
* @param count - count of values in the halo
* @param sendbuf - memory buffer to use to pack values for MPI
* @param ID - 8-bit values on mesh [Nx, Ny, Nz]
*/
void PackID(int *list, int count, signed char *sendbuf, signed char *ID);
/**
* \brief Unpack halo data for 8-bit integer
* @param list - list of values in the halo
* @param count - count of values in the halo
* @param recvbuf - memory buffer containing values recieved by MPI
* @param ID - 8-bit values on mesh [Nx, Ny, Nz]
*/
void UnpackID(int *list, int count, signed char *recvbuf, signed char *ID);
void CommHaloIDs();
//......................................................................................
MPI_Request req1[18], req2[18];
@ -198,6 +268,44 @@ private:
};
template<class TYPE> class PatchData;
enum class DataLocation { CPU, DEVICE };
/**
* \class Patch
*
* @details
* store patch data
*/
class Patch {
public:
//! Empty constructor
Patch() = delete;
//! Copy constructor
Patch( const Patch& ) = delete;
//! Assignment operator
Patch& operator=( const Patch& ) = delete;
//! Return the box for the patch
inline const Box& getBox() const { return d_box; }
//! Create patch data
template<class TYPE>
std::shared_ptr<PatchData<TYPE>> createPatchData( DataLocation location ) const;
private:
Box d_box;
int d_owner;
Domain *d_domain;
};
// Class to hold data on a patch
template<class TYPE>

8
common/ScaLBL.h
View File

@ -234,12 +234,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,double *Vel, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,double *Vel, double *MobilityRatio, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros, double *Perm,double *Vel,double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros, double *Perm,double *Vel, double *MobilityRatio, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);
@ -253,10 +253,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *di
extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np);

103
cpu/GreyscaleColor.cpp
View File

@ -1271,7 +1271,6 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, doubl
-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17);
dist[16*Np+n] = fq;
// q = 17
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)
@ -1341,7 +1340,8 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, doubl
//CP: capillary penalty
// also turn off recoloring for grey nodes
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,
double *Velocity, double *MobilityRatio, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta,
double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
@ -1378,7 +1378,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
/* Corey model parameters */
double Kn_grey,Kw_grey;
double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB;
double GreyDiff; // grey diffusion
double GreyDiff=0.0e-4; // grey diffusion
const double mrt_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802;
@ -1399,7 +1399,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
nB = Den[Np + n];
porosity = Poros[n];
GreyDiff = Perm[n];
//GreyDiff = Perm[n];
perm = 1.0;
W = GreySolidW[n];
Sn_grey = GreySn[n];
@ -1411,6 +1411,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
phi=(nA-nB)/(nA+nB);
// local density
rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
//rho0 *= porosity;
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
tau_eff=tauA_eff + 0.5*(1.0-phi)*(tauB_eff-tauA_eff);
@ -1423,21 +1424,33 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
Krw_grey = 0.0;
if (nA/(nA+nB)<Sn_grey && porosity !=1.0){
perm = Kw_grey;
Krw_grey = Kw_grey;
Swn = 0.0;
}
else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey);
Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero
Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero
Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 4, make sure that W cannot shift to zero
// recompute the effective permeability
perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5));
perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauB-0.5));
//mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5));
}
else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){
perm = Kn_grey;
Krn_grey = Kn_grey;
Swn = 1.0;
}
mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5));
}
/* compute the mobility ratio */
if (porosity != 1.0){
mobility_ratio =(Krn_grey/(tauA-0.5) - Krw_grey/(tauB-0.5))/(Krn_grey/(tauA-0.5) + Krw_grey/(tauB-0.5));
}
else if (phi > 0.0){
mobility_ratio = 1.0;
}
else {
mobility_ratio = -1.0;
}
MobilityRatio[n] = mobility_ratio;
// Get the 1D index based on regular data layout
ijk = Map[n];
@ -2092,7 +2105,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx;
delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx;
jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio);
}
@ -2122,7 +2135,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny;
delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny;
jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio);
}
@ -2153,7 +2166,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz;
delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz;
jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio);
}
@ -2181,7 +2194,8 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//CP: capillary penalty
// also turn off recoloring for grey nodes
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,
double *Velocity, double *MobilityRatio, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
@ -2205,7 +2219,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
/* Corey model parameters */
double Kn_grey,Kw_grey;
double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB;
double GreyDiff; // grey diffusion
double GreyDiff=0.0e-4; // grey diffusion
//double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002)
double porosity;
@ -2235,7 +2249,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
nB = Den[Np + n];
porosity = Poros[n];
GreyDiff = Perm[n];
//GreyDiff = Perm[n];
perm = 1.0;
W = GreySolidW[n];
Sn_grey = GreySn[n];
@ -2248,33 +2262,48 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
// local density
rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
//rho0 *= porosity;
// local relaxation time
tau=tauA + 0.5*(1.0-phi)*(tauB-tauA);
tau_eff=tauA_eff + 0.5*(1.0-phi)*(tauB_eff-tauA_eff);
rlx_setA = 1.f/tau;
rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity
mobility_ratio = 1.0;
Krn_grey = 0.0;
Krw_grey = 0.0;
if (nA/(nA+nB)<Sn_grey && porosity !=1.0){
perm = Kw_grey;
Krw_grey = Kw_grey;
Swn = 0.0;
}
else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey);
Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero
Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero
Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 4, make sure that W cannot shift to zero
// recompute the effective permeability
perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5));
perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauB-0.5));
//mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5));
}
else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){
perm = Kn_grey;
Krn_grey = Kn_grey;
Swn = 1.0;
}
mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5));
}
/* compute the mobility ratio */
if (porosity != 1.0){
mobility_ratio =(Krn_grey/(tauA-0.5) - Krw_grey/(tauB-0.5))/(Krn_grey/(tauA-0.5) + Krw_grey/(tauB-0.5));
}
else if (phi > 0.0){
mobility_ratio = 1.0;
}
else {
mobility_ratio = -1.0;
}
MobilityRatio[n] = mobility_ratio;
// Get the 1D index based on regular data layout
ijk = Map[n];
// COMPUTE THE COLOR GRADIENT
@ -2861,7 +2890,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx;
delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx;
jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio);
}
@ -2887,7 +2916,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny;
delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny;
jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio);
}
@ -2914,7 +2943,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz;
delta = -0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz;
jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio);
}
@ -2964,32 +2993,4 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl
}
}
//extern "C" void ScaLBL_D3Q19_GreyscaleColor_Init(double *dist, double *Porosity, int Np){
// int n;
// double porosity;
// for (n=0; n<Np; n++){
// porosity = Porosity[n];
// if (porosity==0.0) porosity=1.f;
// dist[n] = 0.3333333333333333/porosity;
// dist[Np+n] = 0.055555555555555555/porosity; //double(100*n)+1.f;
// dist[2*Np+n] = 0.055555555555555555/porosity; //double(100*n)+2.f;
// dist[3*Np+n] = 0.055555555555555555/porosity; //double(100*n)+3.f;
// dist[4*Np+n] = 0.055555555555555555/porosity; //double(100*n)+4.f;
// dist[5*Np+n] = 0.055555555555555555/porosity; //double(100*n)+5.f;
// dist[6*Np+n] = 0.055555555555555555/porosity; //double(100*n)+6.f;
// dist[7*Np+n] = 0.0277777777777778/porosity; //double(100*n)+7.f;
// dist[8*Np+n] = 0.0277777777777778/porosity; //double(100*n)+8.f;
// dist[9*Np+n] = 0.0277777777777778/porosity; //double(100*n)+9.f;
// dist[10*Np+n] = 0.0277777777777778/porosity; //double(100*n)+10.f;
// dist[11*Np+n] = 0.0277777777777778/porosity; //double(100*n)+11.f;
// dist[12*Np+n] = 0.0277777777777778/porosity; //double(100*n)+12.f;
// dist[13*Np+n] = 0.0277777777777778/porosity; //double(100*n)+13.f;
// dist[14*Np+n] = 0.0277777777777778/porosity; //double(100*n)+14.f;
// dist[15*Np+n] = 0.0277777777777778/porosity; //double(100*n)+15.f;
// dist[16*Np+n] = 0.0277777777777778/porosity; //double(100*n)+16.f;
// dist[17*Np+n] = 0.0277777777777778/porosity; //double(100*n)+17.f;
// dist[18*Np+n] = 0.0277777777777778/porosity; //double(100*n)+18.f;
// }
//}

34
cpu/Ion.cpp
View File

@ -80,13 +80,14 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
double uEPx,uEPy,uEPz;//electrochemical induced velocity
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
int nr1,nr2,nr3,nr4,nr5,nr6;
@ -124,6 +125,20 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
// q=6
nr6 = neighborList[n+5*Np];
f6 = dist[nr6];
// compute diffusive flux
flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
@ -149,13 +164,14 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
double uEPx,uEPy,uEPz;//electrochemical induced velocity
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
for (n=start; n<finish; n++){
@ -180,6 +196,20 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Veloci
f5 = dist[6*Np+n];
f6 = dist[5*Np+n];
// compute diffusive flux
flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;

4
cpu/Poisson.cpp
View File

@ -132,7 +132,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
f6 = dist[nr6];
Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice speed of sound
Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice squared speed of sound
Ez = (f5-f6)*rlx*4.0;
ElectricField[n+0*Np] = Ex;
ElectricField[n+1*Np] = Ey;
@ -189,7 +189,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_c
Ex = (f1-f2)*rlx*4.0;//NOTE the unit of electric field here is V/lu
Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice speed of sound
Ey = (f3-f4)*rlx*4.0;//factor 4.0 is D3Q7 lattice squared speed of sound
Ez = (f5-f6)*rlx*4.0;
ElectricField[n+0*Np] = Ex;
ElectricField[n+1*Np] = Ey;

1767
cuda/GreyscaleColor.cu
View File

File diff suppressed because it is too large Load Diff

42
cuda/Ion.cu
View File

@ -97,13 +97,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *D
}
}
__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
double uEPx,uEPy,uEPz;//electrochemical induced velocity
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
int nr1,nr2,nr3,nr4,nr5,nr6;
@ -146,6 +147,20 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
nr6 = neighborList[n+5*Np];
f6 = dist[nr6];
// compute diffusive flux
flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
//dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
@ -177,13 +192,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
double uEPx,uEPy,uEPz;//electrochemical induced velocity
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
int S = Np/NBLOCKS/NTHREADS + 1;
@ -212,6 +228,20 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *V
f5 = dist[6*Np+n];
f6 = dist[5*Np+n];
// compute diffusive flux
flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
//dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
@ -330,10 +360,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAodd_Ion<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Den,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
dvc_ScaLBL_D3Q7_AAodd_Ion<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
@ -342,10 +372,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAeven_Ion<<<NBLOCKS,NTHREADS >>>(dist,Den,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
dvc_ScaLBL_D3Q7_AAeven_Ion<<<NBLOCKS,NTHREADS >>>(dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){

BIN
docs/source/_static/images/AA-stream.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 332 KiB

BIN
docs/source/_static/images/data-structures-v0.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 165 KiB

BIN
docs/source/_static/images/domain-decomp.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 272 KiB

BIN
docs/source/_static/images/lbpm-visit-workflow-i.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 516 KiB

BIN
docs/source/_static/images/lbpm-visit-workflow-ii.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 371 KiB

BIN
docs/source/_static/images/lbpm-visit-workflow-iii.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 455 KiB

BIN
docs/source/_static/images/lbpm-visit-workflow-iv.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 1.9 MiB

119
docs/source/developerGuide/designOverview/dataStructures.rst
View File

@ -3,4 +3,121 @@ Data Structures
===============
LBPM includes a variety of generalized data structures to facilitate the implementation
of different lattice Boltzmann models.
of different lattice Boltzmann models. These core data structures are designed so that
they are agnostic to the particular physics of a model. Many different lattice Boltzmann
schemes can be constructed using the same basic data structures. The core data structures
are primarily designed to
1) prepare 3D image data for distributed memory simulation based on MPI
2) develop common data structures to facilitate the local streaming step
3) generate and store information needed to exchange data between MPI processes
4) develop common routines for parallel communication for lattice Boltzmann methods.
By using the core data structures, it is possible to develop new physical models
without having to revisit the data structures or parallel algorithms on which a model is
built. Understanding and using the core data structures is very useful if you intend
to build new physical models within LBPM. In most cases, only the collision step
will need to be developed to implement a particular physical model.
------------------
Domain
------------------
The ``Domain`` data structure includes information needed to carry out generic
parallel computations on 3D image data. LBPM is designed to support parallel
simulation where the image data is distributed over a potentially large number of
processors. Each processor recieves an equal size chunk of the input image,
which is specified based on the input database file in the form ``n = nx, ny, nz``.
The total image size is specified as ``N = Nx, Ny, Nz``, which is needed to read
the input image. In the typical case, the input image data will be read only by MPI
rank 0, and sub-domains will be assigned and distributed to each MPI process.
The regular process grid is specified as ``nproc = Px, Py, Pz``. If the entire image
is to be distributed across the processors, the values must be chosen such that
``Nx = nx*Px``, ``Ny = ny*Py`` and ``Nz = nz*Pz``. The basic purpose for the
data structures contained in the ``Domain`` class is to retain basic information
about the domain structure to support parallel computations in MPI. Code within the
``analysis/`` directory will generally rely on these data structures to determine
how to perform computations in parallel. For GPU-based application, information
contained within the ``Domain`` structure is usually retained only the CPU.
Lattice Boltzmann algorithms rely on dedicated data structures within the ``ScaLBL``
class to perform MPI computations, which are designed for these applications
specifically.
.. figure:: ../../_static/images/domain-decomp.png
:width: 600
:alt: domain decomposition
Domain decomposition used by LBPM to distribute 3D image data across processors.
Each processor recieves an equal size block of the input image. The Domain class
organizes MPI communications with adjacent blocks based on the processor layout.
------------------
ScaLBL
------------------
While the ``Domain`` data structures retain essential information needed to
carry out generic types of computations based on regular 3D image files,
``ScaLBL`` data structures are designed specifically for lattice Boltzmann methods.
The main purposes for these data structures are to:
1) allow for overlap between computations that do not depend on MPI communication and the
MPI communications required by LBMs (e.g. to carry out the streaming step, perform
halo exchanges, etc.)
2) reduce the computational and memory requirements by excluding image sites where flow
does not occur; and
3) to facilitate efficient data access patterns based on the memory layout used by
data structures.
These are critical steps to provide good parallel performance and make efficient use
of GPUs.
In many cases flow will occur only in a fraction of the voxels from the original 3D input
image. Since LBM schemes store a lot of data for each active voxels (e.g. 19 double precision
values for a single D3Q19 model), significant memory savings can be realized by excluding
immobile voxels from the data structure. Furthermore, it is useful to distiguish between
the interior lattice sites (which do not require MPI communications to update) and the exterior
lattice sites (which must wait for MPI communications to complete before the local update can
be completed). Within LBPM it is almost always advantageous to overlap interior computations with
MPI communications, so that the latter costs can be hidden behind the computational density
of the main lattice Boltzmann kernels. As the local sub-domain size increases, the number of
interior sites will tend to be significantly larger than the number of exterior sites.
The data layout is summarized in the image below.
.. figure:: ../../_static/images/data-structures-v0.png
:width: 600
:alt: ScaLBL data structure layout
Data structure used by ScaLBL to support lattice Boltzmann methods.
Regions of the input image that are occupied by solid are excluded.
Interior and exterior lattice sites are stored separately so that
computations can efficiently overlap with MPI communication.
Another factor in the data structure design are the underlying algorithms used by the LBM
to carry out the streaming step. There are at least a half dozen distinct algorithms that
can be used to perform streaming for lattice Boltzmann methods, each with their own advantages
and disadvantages. Generally, the choice of streaming algorithm should reduce the overall memory
footprint while also avoiding unnecessary memory accesses. LBPM uses the AA algorithm so
that a single array is needed to perform the streaming step. The AA algorithm relies on the symmetry
of the discrete velocity set used by LBMs to implement an efficient algorithm. Conceptually,
we can easily see that for each distribution that needs to be accessed to perform streaming,
there is an opposite distribution. The AA algorithm proceeds by identifying these distributions
and exchanging their storage locations in memory, as depicted in the image below. The data
access pattern will therefore alternate between even and odd timesteps. Within LBPM, a ``NeighborList``
is constructed to store the memory location for the accompanying pair for each distribution.
If the neighboring site is in the solid, the neighbor is specified such that the halfway
bounceback rule will applied automatically. In this way, the streaming step can be performed
very efficiently on either GPU or CPU.
.. figure:: ../../_static/images/AA-stream.png
:width: 600
:alt: AA streaming pattern
Data access pattern for AA streaming algorithm used by LBPM. The location to read
and write distributions alternates between even and odd timesteps. A solid bounceback
rule is built into the data structure.

7
docs/source/developerGuide/doxygen/Domain.rst Normal file
View File

@ -0,0 +1,7 @@
============================================
Domain
============================================
.. doxygenfile:: Domain.h
:project: LBPM Doxygen

2
docs/source/developerGuide/index.rst
View File

@ -17,6 +17,8 @@ into the framework.
doxygen/ScaLBL.rst
doxygen/Domain.rst
doxygen/Models.rst
doxygen/Analysis.rst

5
docs/source/examples/color/steadyState.rst Normal file
View File

@ -0,0 +1,5 @@
*******************
Steady-state flow
*******************
In this example we simulate a steady-state flow with a constant driving force.

52
docs/source/examples/domain/domain.rst Normal file
View File

@ -0,0 +1,52 @@
*****************
Input Domain
*****************
LBPM provides a flexible framework to ingest 3D image data.
To illustrate the basic capabilities, this tutorial considers a quasi-2D
flow cell. Source files for the example are included in the LBPM repository
in the directory ``examples/DiscPack``. A simple python code is included
to set up the flow domain.
Based on LBPM convention, external boundary conditions are applied in the
z-direction. This means that the domain should be set up so that the direction
you want to set boundary conditions is aligned with the z-axis. For the quasi-2D
example, a depth of ``3`` voxels is used for the x-direction. *Based on LBPM
internal data structures at least three voxels must be provided in each direction*
The specified domain decomposition must also observe this rule.
Image data is stored internally within LBPM as signed 8-bit binary data. This means that
up to 256 labels can be provided in the input image. LBPM convention takes all
non-positive labels to be immobile (treated as solid). In this example, the solid regions
are assigned a value of ``0``. It is possible to provide up to ``128`` different labels
for the solid. Also, note that python supports only the unsigned 8-bit datatype. For the unsigned data
type, labels assigned values ``128,...255`` in python will correspond to labels
``-127,...-1`` when read in as type ``signed char`` within LBPM.
.. code:: python
import numpy as np
import matplotlib.pylab as plt
import pandas as pd
# Set the size of the domain
Nx=3
Ny=128
Nz=128
D=pd.read_csv("discs.csv",sep=" ")
ID = np.ones(Nx*Ny*Nz,dtype='uint8')
ID.shape = (Nz,Ny,Nx)
# Set the solid labels
for idx in range(len(D)):
cx=D['cx'][idx] / dx
cy=D['cy'][idx] /dx
r=D['r'][idx] /dx
for i in range(0,Nz):
for j in range(0,Ny):
if ( (cx-i)*(cx-i) + (cy-j)*(cy-j) < r*r ):
ID[i,j,0] = 0
ID[i,j,1] = 0
ID[i,j,2] = 0
# write input file to disc
ID.tofile("discs_3x128x128.raw")

106
docs/source/examples/morphology/morphOpen.rst Normal file
View File

@ -0,0 +1,106 @@
*****************************
Morphological pre-processors
*****************************
LBPM includes morphological pre-processing tools as utility functions.
It is often useful to generate initial conditions for a 2-phase flow simulation based on a morphological approach. In particular, morphological tools can be used to provide a physical reasonable initial condition in cases where direct experimental observations are not available. These initial configurations are compatible with any of the 2-phase simulation protocols used by lbpm_color_simulator. These initialization approaches alter the fluid labels within the input files, writing a new file with the new morphologically assigned labels.
There are two main morphological pre-processors in LBPM
* ``lbpm_morphdrain_pp`` -- initialize fluid configuration based on morphological drainage
* ``lbpm_morphopen_pp`` -- initialize fluid configuration based on morphological opening
Here we demonstrate ``lbpm_morphdrain_pp`` because it offers the most information. Although it is not perfect, the morphological drainage operation does a good job of approximating configurations observed along primary drainage processes performed under water-wet conditions. A limitation is that fluid trapped in the corners will not stop the morphological operation from eroding it. This should not discourage you too much -- morphological tools are very practical and can save you a lot of time! It is also a good thing to be skeptical.
Since the morphological operation works on the input domain, associated parameters are added to the ``Domain`` section of the input file. Here we will set a target saturation Sw = 0.20, which will run the morphological drainage operation until the fluid labeled as 2 occupies 20% of the pore space or less. For the case considered in
``example/DiscPack`` we specify the following information in the input file
.. code:: c
Domain {
Filename = "discs_3x128x128.raw"
ReadType = "16bit" // data type
N = 3, 128, 128 // size of original image
nproc = 1, 2, 2 // process grid
n = 3, 64, 64 // sub-domain size
voxel_length = 7.0 // voxel length (in microns)
ReadValues = 0, 1, 2 // labels within the original image
WriteValues = 0, 2, 2 // associated labels to be used by LBPM
InletLayers = 0, 0, 6 // specify 6 layers along the z-inlet
BC = 0 // fully periodic BC
Sw = 0.20 // target saturation for morphological tools
}
Once this has been set, we launch lbpm_morphdrain_pp in the same way as other parallel tools
.. code:: bash
mpirun -np 4 $LBPM_BIN/lbpm_morphdrain_pp input.db
Successful output looks like the following
.. code:: bash
Performing morphological opening with target saturation 0.500000
voxel length = 1.000000 micron
voxel length = 1.000000 micron
Input media: discs_3x128x128.raw
Relabeling 3 values
oldvalue=0, newvalue =0
oldvalue=1, newvalue =2
oldvalue=2, newvalue =2
Dimensions of segmented image: 3 x 128 x 128
Reading 8-bit input data
Read segmented data from discs_3x128x128.raw
Label=0, Count=11862
Label=1, Count=37290
Label=2, Count=0
Distributing subdomains across 4 processors
Process grid: 1 x 2 x 2
Subdomain size: 3 x 64 x 64
Size of transition region: 0
Media porosity = 0.758667
Initialized solid phase -- Converting to Signed Distance function
Volume fraction for morphological opening: 0.758667
Maximum pore size: 116.773801
1.000000 110.935111
1.000000 105.388355
1.000000 100.118937
1.000000 95.112990
1.000000 90.357341
1.000000 85.839474
1.000000 81.547500
1.000000 77.470125
1.000000 73.596619
1.000000 69.916788
1.000000 66.420949
1.000000 63.099901
1.000000 59.944906
1.000000 56.947661
1.000000 54.100278
1.000000 51.395264
1.000000 48.825501
1.000000 46.384226
1.000000 44.065014
1.000000 41.861764
1.000000 39.768675
1.000000 37.780242
1.000000 35.891230
1.000000 34.096668
1.000000 32.391835
0.805149 30.772243
0.727353 29.233631
0.719791 27.771949
0.714883 26.383352
0.710861 25.064184
0.637249 23.810975
0.444570 22.620426
Final void fraction =0.444570
Final critical radius=22.620426
Writing ID file
Writing file to: discs_3x128x128.raw.morphdrain.raw
The final configuration can be visualized in python by loading the output file
``discs_3x128x128.raw.morphdrain.raw``.

15
docs/source/examples/running.rst
View File

@ -1,10 +1,21 @@
==============
LBPM examples
Examples
==============
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.
understanding the LBPM input file structure. The examples included provide
a basic introduction to working with LBPM.
.. toctree::
:glob:
:maxdepth: 2
domain/*
morphology/*
color/*

2
docs/source/index.rst
View File

@ -12,9 +12,9 @@ LBPM -- Documentation
:caption: Contents:
install
examples/*
userGuide/*
developerGuide/*
examples/*
publications/*
Indices and tables

6
docs/source/userGuide/IO/fileformat.rst
View File

@ -5,9 +5,9 @@ I/O conventions for LBPM
There are three main kinds of output file that are supported by LBPM.
* CSV files --
* CSV files -- space-delimited CSV files are used by the internal analysis framework
* formatted binary files --
* formatted binary files -- SILO and HDF5 formats are supported for visualization data
* unformatted binary files --
* unformatted binary files -- ``.raw`` extension

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

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

2
docs/source/userGuide/index.rst
View File

@ -10,8 +10,6 @@ Welcome to the LBPM user guide.
models/*
analysis/*
visualization/*
IO/*

5
docs/source/userGuide/models/color/index.rst
View File

@ -266,3 +266,8 @@ the inlet or outlet, the ``Domain`` section of the database may specify the foll
- ``InletLayerPhase = 2`` -- establish a reservoir of component B at the inlet
- ``OutletLayerPhase = 1`` -- establish a reservoir of component A at the outlet
****************
Example data
****************
Example data can be downloaded from https://www.digitalrocksportal.org/projects/326

76
docs/source/userGuide/models/color/protocols/centrifuge.rst
View File

@ -30,14 +30,76 @@ This is because a positive body force will favor a larger saturation of fluid A
(positive capillary pressure ) whereas a negative body force will favor a lower
saturation of fluid A (negative capillary pressure).
The simplest way to infer the capillary pressure is based on consideration of the average
fluid pressures, which are logged to the output files ``timelog.csv`` and ``subphase.csv``.
In the units of the lattice Boltzmann simulation, the interfacial tension is given
as :math:`\gamma_{wn} = 6 \alpha`. Suppose that the physical interfacial tension is given by
:math:`\gamma_{wn}^\prime`, provided in units of Pa-m. The capillary pressure in pascal will
then be given by
To enable the ``centrifuge`` protocol such that the pressure of fluid B is higher than
fluid A, the following keys can be set. Increasing the body force will lead to a larger
capillary pressure
.. math::
:nowrap:
$$
p_c^\prime = \frac{\gamma_{wn}^\prime (p_n - p_w)}{\gamma_{wn} \Delta x}
$$
where :math:`\Delta x` is the voxel length in meters.
To enable the ``centrifuge`` protocol such that the effective pressure of fluid B is higher
than fluid A, the input file can be specified as below. Increasing the body force will lead to
a larger capillary pressure.
.. code-block:: c
Color {
protocol = "centrifuge"
timestepMax = 1000000 // maximum timtestep
alpha = 0.005 // controls interfacial tension
rhoA = 1.0 // controls the density of fluid A
rhoB = 1.0 // controls the density of fluid B
tauA = 0.7 // controls the viscosity of fluid A
tauB = 0.7 // controls the viscosity of fluid B
F = 0, 0, -1.0e-5 // body force
din = 1.0 // inlet density (controls pressure)
dout = 1.0 // outlet density (controls pressure)
WettingConvention = "SCAL" // convention for sign of wetting affinity
ComponentLabels = 0, -1, -2 // image labels for solid voxels
ComponentAffinity = 1.0, 1.0, 0.6 // controls the wetting affinity for each label
Restart = false
}
Domain {
Filename = "Bentheimer_LB_sim_intermediate_oil_wet_Sw_0p37.raw"
ReadType = "8bit" // data type
N = 900, 900, 1600 // size of original image
nproc = 2, 2, 2 // process grid
n = 200, 200, 200 // sub-domain size
offset = 300, 300, 300 // offset to read sub-domain
voxel_length = 1.66 // voxel length (in microns)
ReadValues = -2, -1, 0, 1, 2 // labels within the original image
WriteValues = -2, -1, 0, 1, 2 // associated labels to be used by LBPM
BC = 3 // boundary condition type (0 for periodic)
}
Analysis {
analysis_interval = 1000 // logging interval for timelog.csv
subphase_analysis_interval = 5000 // loggging interval for subphase.csv
visualization_interval = 100000 // interval to write visualization files
N_threads = 4 // number of analysis threads (GPU version only)
restart_interval = 1000000 // interval to write restart file
restart_file = "Restart" // base name of restart file
}
Visualization {
write_silo = true // write SILO databases with assigned variables
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
save_phase_field = true // save phase field within SILO database
save_pressure = false // save pressure field within SILO database
save_velocity = false // save velocity field within SILO database
}
FlowAdaptor {
}
.. code-block:: bash
protocol = "centrifuge"
F = 0, 0, -1.0e-5

48
docs/source/userGuide/models/color/protocols/coreFlooding.rst
View File

@ -33,7 +33,7 @@ The volumetric flow rate is related to the fluid velocity based on
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
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
@ -59,5 +59,51 @@ will necessarily restrict the LB parameters that can be used, which are ultimate
manifested as a constraint on the size of a timestep.
.. code-block:: c
Color {
protocol = "core flooding"
capillary_number = 1e-4 // capillary number for the displacement
timestepMax = 1000000 // maximum timtestep
alpha = 0.005 // controls interfacial tension
rhoA = 1.0 // controls the density of fluid A
rhoB = 1.0 // controls the density of fluid B
tauA = 0.7 // controls the viscosity of fluid A
tauB = 0.7 // controls the viscosity of fluid B
F = 0, 0, 0 // body force
WettingConvention = "SCAL" // convention for sign of wetting affinity
ComponentLabels = 0, -1, -2 // image labels for solid voxels
ComponentAffinity = 1.0, 1.0, 0.6 // controls the wetting affinity for each label
Restart = false
}
Domain {
Filename = "Bentheimer_LB_sim_intermediate_oil_wet_Sw_0p37.raw"
ReadType = "8bit" // data type
N = 900, 900, 1600 // size of original image
nproc = 2, 2, 2 // process grid
n = 200, 200, 200 // sub-domain size
offset = 300, 300, 300 // offset to read sub-domain
voxel_length = 1.66 // voxel length (in microns)
ReadValues = -2, -1, 0, 1, 2 // labels within the original image
WriteValues = -2, -1, 0, 1, 2 // associated labels to be used by LBPM
BC = 4 // boundary condition type (0 for periodic)
}
Analysis {
analysis_interval = 1000 // logging interval for timelog.csv
subphase_analysis_interval = 5000 // loggging interval for subphase.csv
visualization_interval = 100000 // interval to write visualization files
N_threads = 4 // number of analysis threads (GPU version only)
restart_interval = 1000000 // interval to write restart file
restart_file = "Restart" // base name of restart file
}
Visualization {
write_silo = true // write SILO databases with assigned variables
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
save_phase_field = true // save phase field within SILO database
save_pressure = false // save pressure field within SILO database
save_velocity = false // save velocity field within SILO database
}
FlowAdaptor {
}

108
docs/source/userGuide/models/color/protocols/fractionalFlow.rst
View File

@ -6,12 +6,110 @@ 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.
the local rate of flow is high.
.. code-block:: bash
Enabling the fractional flow protocol will automatically set the following default settings:
protocol = "fractional flow"
* ``fractional_flow_increment = 0.05`` - target change in saturation between steady points
* ``endpoint_threshold = 0.1`` - termination criterion based on the relative flow rates of fluids
* ``mass_fraction_factor = 0.006`` - percentage change to the mass fraction (per iteration)
* ``fractional_flow_epsilon = 5.0e-6`` - control the threshold velocity :math:`U_\epsilon`
* ``skip_timesteps = 50000`` - timesteps to spend in adaptive part of algorithm
* ``min_steady_timesteps = 1000000`` - minimum number of timesteps per steady point
* ``max_steady_timesteps = 1000000`` - maximum number of timesteps per steady point
Any of these values can be manually overriden by setting the appropriate key to a custom value
within the ``FlowAdaptor`` section of the input file database. The fractional flow
will enforce periodic boundary conditions ``BC = 0``, and is not compatible with other
boundary condition routines. A warning message will be printed if ``BC`` is set to a value
other than ``0``.
The basic idea for the fractional flow algorithm is to define an algorithm to modify the
fluid saturation that will:
(1) minimize the introduction of end effects and gradients in saturation
(2) respect the structure and connectivity of the fluid phases to the maximum extent
(3) minimize the transient disruption to the flow behavior so that fewer timesteps are required to reach steady
state
The strategy used to accomplish these objectives is to change saturation by either increasing
or decreasing the mass within each fluid, depending on the sign of ``mass_fraction_factor``.
If ``mass_fraction_factor`` is positive, the algorithm will remove mass from fluid A and add
mass to fluid B. The opposite will happen if ``mass_fraction_factor`` is negative. The algorithm
will selectively add mass to voxels where the local flow rate is high, since these will generally
correspond to the mobile (i.e. not disconnected) fluid regions. A local weighting function is determined
from the velocity :math:`\mathbf{u}_i`
.. math::
:nowrap:
$$
W_i = \frac{ U_\epsilon + |\mathbf{u}_i|}{U_\epsilon + \max{|\mathbf{u}_i|}}
$$
where :math:`\max{|\mathbf{u}_i|}` is the maximum flow speed within fluid :math:`i` and
:math:`U_\epsilon` is a threshold speed that is set to minimize the influence of spurious
currents on the mass seeding algorithm. The sum of the weighting function is used to normalize
the local weights so that the added mass will match the value specified by ``mass_fraction_factor``.
If the flow is slower than :math:`U_\epsilon`, the algorithm will tend to add mass evenly to the system.
For example, if the water is only present in films that flow very slowly, then mass will
be evenly seeded throughout entire water film. Alternatively, if one or both fluids
flows through distinct channels, the mass will be disproportionately added to these
channels where the rate of flow is high. As system relaxes, the mass will redistribute
spatially, causing a change to the fluid saturation.
.. code-block:: c
Color {
protocol = "fractional flow"
capillary_number = 1e-4 // capillary number for the displacement
timestepMax = 1000000 // maximum timtestep
alpha = 0.005 // controls interfacial tension
rhoA = 1.0 // controls the density of fluid A
rhoB = 1.0 // controls the density of fluid B
tauA = 0.7 // controls the viscosity of fluid A
tauB = 0.7 // controls the viscosity of fluid B
F = 0, 0, 1.0e-5 // body force
WettingConvention = "SCAL" // convention for sign of wetting affinity
ComponentLabels = 0, -1, -2 // image labels for solid voxels
ComponentAffinity = 1.0, 1.0, 0.6 // controls the wetting affinity for each label
Restart = false
}
Domain {
Filename = "Bentheimer_LB_RelPerm_intermediate_oil_wet_Sw_0p37.raw"
ReadType = "8bit" // data type
N = 900, 900, 1600 // size of original image
nproc = 2, 2, 2 // process grid
n = 200, 200, 200 // sub-domain size
offset = 300, 300, 300 // offset to read sub-domain
InletLayers = 0, 0, 6 // number of mixing layers at the inlet
OutletLayers = 0, 0, 6 // number of mixing layers at the outlet
voxel_length = 1.66 // voxel length (in microns)
ReadValues = -2, -1, 0, 1, 2 // labels within the original image
WriteValues = -2, -1, 0, 1, 2 // associated labels to be used by LBPM
BC = 0 // boundary condition type (0 for periodic)
}
Analysis {
analysis_interval = 1000 // logging interval for timelog.csv
subphase_analysis_interval = 5000 // loggging interval for subphase.csv
visualization_interval = 100000 // interval to write visualization files
N_threads = 4 // number of analysis threads (GPU version only)
restart_interval = 1000000 // interval to write restart file
restart_file = "Restart" // base name of restart file
}
Visualization {
write_silo = true // write SILO databases with assigned variables
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
save_phase_field = true // save phase field within SILO database
save_pressure = false // save pressure field within SILO database
save_velocity = false // save velocity field within SILO database
}
FlowAdaptor {
min_steady_timesteps = 100000 // minimum number of timesteps per steady point
max_steady_timesteps = 250000 // maximum number of timesteps per steady point
mass_fraction_factor = 0.006 // controls the rate of mass seeding in adaptive step
fractional_flow_increment = 0.05 // saturation change after each steady point
endpoint_threshold = 0.1 // endpoint exit criterion (based on flow rates)
}

52
docs/source/userGuide/models/color/protocols/imageSequence.rst
View File

@ -21,7 +21,53 @@ To enable the image sequence protocol, the following keys should be set
within the ``Color`` section of the input database
.. code-block:: bash
Color {
protocol = "image sequence"
image_sequence = "Bentheimer_LB_RelPerm_intermediate_oil_wet_Sw_0p37.raw", "Bentheimer_LB_RelPerm_intermediate_oil_wet_Sw_0p72.raw"
// comma separated list of image names
capillary_number = 1e-4 // capillary number for the displacement
timestepMax = 1000000 // maximum timtestep
alpha = 0.005 // controls interfacial tension
rhoA = 1.0 // controls the density of fluid A
rhoB = 1.0 // controls the density of fluid B
tauA = 0.7 // controls the viscosity of fluid A
tauB = 0.7 // controls the viscosity of fluid B
F = 0, 0, 0 // body force
WettingConvention = "SCAL" // convention for sign of wetting affinity
ComponentLabels = 0, -1, -2 // image labels for solid voxels
ComponentAffinity = 1.0, 1.0, 0.6 // controls the wetting affinity for each label
Restart = false
}
Domain {
Filename = "Bentheimer_LB_RelPerm_intermediate_oil_wet_Sw_0p37.raw"
ReadType = "8bit" // data type
N = 900, 900, 1600 // size of original image
nproc = 2, 2, 2 // process grid
n = 200, 200, 200 // sub-domain size
offset = 300, 300, 300 // offset to read sub-domain
voxel_length = 1.66 // voxel length (in microns)
ReadValues = -2, -1, 0, 1, 2 // labels within the original image
WriteValues = -2, -1, 0, 1, 2 // associated labels to be used by LBPM
BC = 0 // boundary condition type (0 for periodic)
}
Analysis {
analysis_interval = 1000 // logging interval for timelog.csv
subphase_analysis_interval = 5000 // loggging interval for subphase.csv
visualization_interval = 100000 // interval to write visualization files
N_threads = 4 // number of analysis threads (GPU version only)
restart_interval = 1000000 // interval to write restart file
restart_file = "Restart" // base name of restart file
}
Visualization {
write_silo = true // write SILO databases with assigned variables
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
save_phase_field = true // save phase field within SILO database
save_pressure = false // save pressure field within SILO database
save_velocity = false // save velocity field within SILO database
}
FlowAdaptor {
min_steady_timesteps = 100000 // minimum number of timesteps per steady point
max_steady_timesteps = 250000 // maximum number of timesteps per steady point
}
protocol = "image sequence"
image_sequence = "image1.raw", "image2.raw"

50
docs/source/userGuide/models/color/protocols/shellAggregation.rst
View File

@ -16,3 +16,53 @@ protocol is described by Wang et al. ( https://doi.org/10.1016/j.jcp.2019.108966
.. code-block:: c
Color {
protocol = "shell aggregation"
capillary_number = 1e-4 // capillary number for the displacement
timestepMax = 1000000 // maximum timtestep
alpha = 0.005 // controls interfacial tension
rhoA = 1.0 // controls the density of fluid A
rhoB = 1.0 // controls the density of fluid B
tauA = 0.7 // controls the viscosity of fluid A
tauB = 0.7 // controls the viscosity of fluid B
F = 0, 0, 0 // body force
WettingConvention = "SCAL" // convention for sign of wetting affinity
ComponentLabels = 0, -1, -2 // image labels for solid voxels
ComponentAffinity = 1.0, 1.0, 0.6 // controls the wetting affinity for each label
Restart = false
}
Domain {
Filename = "Bentheimer_LB_sim_intermediate_oil_wet_Sw_0p37.raw"
ReadType = "8bit" // data type
N = 900, 900, 1600 // size of original image
nproc = 2, 2, 2 // process grid
n = 200, 200, 200 // sub-domain size
offset = 300, 300, 300 // offset to read sub-domain
InletLayers = 0, 0, 6 // number of mixing layers at the inlet
OutletLayers = 0, 0, 6 // number of mixing layers at the outlet
voxel_length = 1.66 // voxel length (in microns)
ReadValues = -2, -1, 0, 1, 2 // labels within the original image
WriteValues = -2, -1, 0, 1, 2 // associated labels to be used by LBPM
BC = 0 // boundary condition type (0 for periodic)
}
Analysis {
analysis_interval = 1000 // logging interval for timelog.csv
subphase_analysis_interval = 5000 // loggging interval for subphase.csv
visualization_interval = 100000 // interval to write visualization files
N_threads = 4 // number of analysis threads (GPU version only)
restart_interval = 1000000 // interval to write restart file
restart_file = "Restart" // base name of restart file
}
Visualization {
write_silo = true // write SILO databases with assigned variables
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
save_phase_field = true // save phase field within SILO database
save_pressure = false // save pressure field within SILO database
save_velocity = false // save velocity field within SILO database
}
FlowAdaptor {
}

87
docs/source/userGuide/models/greyscale/greyscale.rst
View File

@ -1,7 +1,90 @@
=============================================
Greyscale model model
Greyscale 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.
solution in open regions. To use the greyscale model, the input image should be segmented
into "grey" and open regions. Each particular "grey" label should be assigned both
a porosity and permeability value.
A typical command to launch the LBPM color simulator is as follows
```
mpirun -np $NUMPROCS lbpm_greyscale_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.
***************************
Model parameters
***************************
The essential model parameters for the single phase greyscale model are
- ``tau`` -- control the fluid viscosity -- :math:`0.7 < \tau < 1.5`
The kinematic viscosity is given by
***************************
Model formulation
***************************
A D3Q19 LBE is constructed to describe the momentum transport
.. math::
:nowrap:
$$
f_q(\bm{x}_i + \bm{\xi}_q \delta t,t + \delta t) - f_q(\bm{x}_i,t) =
\sum^{Q-1}_{k=0} M^{-1}_{qk} S_{kk} (m_k^{eq}-m_k) + \sum^{Q-1}_{k=0} M^{-1}_{qk} (1-\frac{S_{kk}}{2}) \hat{F}_q\;,
$$
The force is imposed based on the construction developed by Guo et al
.. math::
:nowrap:
$$
F_i = \rho_0 \omega_i \left[\frac{\bm{e}_i \cdot \bm{a}}{c_s^2} +
\frac{\bm{u} \bm{a}:(\bm{e}_i \bm{e}_i -c_s^2 \mathcal{I})}{\epsilon c_s^4} \right] ,
$$
where :math:`c_s^2 = 1/3` is the speed of sound for the D3Q19 lattice.
The acceleration includes contributions due to the external driving force :math:`\bm{g}`
as well as a drag force due to the permeability :math:`K` and flow rate :math:`\bm{u}` with the
porosity :math:`\epsilon` and viscosity :math:`\nu` determining the net forces acting within
a grey voxel
.. math::
:nowrap:
$$
\bm{a} = - \frac{\epsilon \nu}{K} \bm{u} + \bm{F}_{cp}/\rho_0 + \epsilon \bm{g},
$$
The flow velocity is defined as
.. math::
:nowrap:
$$
\rho_0 \bm{u} = \sum_i \bm{e}_i f_i + \frac{\delta t}{2} \rho_0 \bm{a}.
$$
Combining the previous expressions,
.. math::
:nowrap:
$$
\bm{u} = \frac{\frac{1}{\rho_0}\sum_i \bm{e}_i f_i + \frac{\delta t}{2}\epsilon \bm{g} +
\frac{\delta t}{2} \frac{\bm{F}_{cp}}{\rho_0}}{1+ \frac{\delta t}{2} \frac{\epsilon \nu}{K}}
$$

330
docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst Normal file
View File

@ -0,0 +1,330 @@
=============================================
Greyscale Color Model
=============================================
The LBPM greyscale lattice Boltzmann model is constructed to approximate the
solution of the Darcy-Brinkman equations in grey regions coupled to a color model implementation
solution in open regions. A simple constitutive form is used to model the relative
permeability in the grey regions,
***************************
Parameters
***************************
The essential model parameters for the color model are
- ``alpha`` -- control the interfacial tension between fluids -- :math:`0 < \alpha < 0.01`
- ``beta`` -- control the width of the interface -- :math:`\beta < 1`
- ``tauA`` -- control the viscosity of fluid A -- :math:`0.7 < \tau_A < 1.5`
- ``tauB`` -- control the viscosity of fluid B -- :math:`0.7 < \tau_B < 1.5`
- ``rhoA`` -- control the viscosity of fluid A -- :math:`0.05 < \rho_A < 1.0`
- ``rhoB`` -- control the viscosity of fluid B -- :math:`0.05 < \rho_B < 1.0`
- ``greyscale_endpoint_A`` -- control the endpoint for greyscale components -- :math:`S_b^{r,a}`
- ``greyscale_endpoint_B`` -- control the endpoint for greyscale components -- :math:`S_b^{r,b}`
****************************
Formulation
****************************
The greyscale color model extends the conventional two-fluid color model such that flow through
micro-porous "grey" voxels can also be treated. Extensions to the formulation are made to:
(1) implement momentum transport behavior consistent with the single phase greyscale model
(2) adapt the re-coloring term to allow for transport through microporosity.
Mass transport LBEs are constructed to model the behavior for each fluid. Within the open pore
regions, mass transport LBEs recover the behavior of the standard two-fluid color model.
Within grey voxels the standard recoloring term is disabled so that the mass transport behavior
can be described by a different rule within the microporous regions. The endpoints are
specified based on the saturation of fluid B at which fluid A ceases to percolate, :math:`S_b^{r,a}`
and the saturation of fluid B at which fluid B ceases to percolate, :math:`S_b^{r,b}`.
The endpoints should be independently specified for each class of microporosity that is labeled
in the input image. Between the endpoint values, the effective permeability is paramaterized based on the value
.. math::
:nowrap:
$$
S_{ab} = \frac{S_b - S_b^{r,b}}{S_b^{r,a} - S_b^{r,b}}
$$
At the endpoints, the effective permeability is provided as an input parameter for
each fluid. When :math:`S_b=S_b^{r,b}` the effective permeability of fluid B is
zero and :math:`K_a=K^{r,a}`. Between the endpoints the Corey model predicts the
effective permeability for fluid A according to
.. math::
:nowrap:
$$
K_a = K^{r,a} (1-S_{ab})^{\lambda^a}
$$
where :math:`\lambda^a=2` is the Corey exponent. Likewise,
When :math:`S_b=S_b^{r,a}` the effective permeability for fluid A will be zero,
and :math:`K_b=K^{r,b}` with
.. math::
:nowrap:
$$
K_b = K^{r,b} S_{ab}^{\lambda^b}
$$
Two LBEs are constructed to model the mass transport,
.. math::
:nowrap:
$$
A_q(\bm{x} + \bm{\xi}_q \delta t, t+\delta t) = w_q N_a \Big[1 + \frac{\bm{u} \cdot \bm{\xi}_q}{c_s^2}
+ \beta \frac{N_b}{N_a+N_b} \bm{n} \cdot \bm{\xi}_q\Big] \;
$$
.. math::
:nowrap:
$$
B_q(\bm{x} + \bm{\xi}_q \delta t, t+\delta t) =
w_q N_b \Big[1 + \frac{\bm{u} \cdot \bm{\xi}_q}{c_s^2}
- \beta \frac{N_a}{N_a+N_b} \bm{n} \cdot \bm{\xi}_q\Big]\;,
$$
The number density for each fluid is obtained from the sum of the mass transport distributions
.. math::
:nowrap:
$$
N_a = \sum_q A_q\;, \quad N_b = \sum_q B_q\;
$$
The phase indicator field is then defined as
.. math::
:nowrap:
$$
\phi = \frac{N_a-N_b}{N_a+N_b}
$$
The recoloring step incorporates the standard color model
rule to model anti-diffusion at the interface within the open pores. Within grey regions,
the anti-diffusion term must be disabled, since within grey voxels the length scale for fluid features
is smaller than the interface width produced from the color model. Within grey voxels the
two fluids are permitted to freely mix between the endpoints. Beyond the endpoints, the recoloring
term is used to drive spontaneous imbibition into the grey voxels
..math::
:nowrap:
$$
R_c =
$$
The fluid density and kinematic viscosity are determined based on linear interpolation
.. math::
:nowrap:
$$
\rho_0 = \frac{(1+\phi) \rho_n}{2}+ \frac{(1-\phi) \rho_w}{2} \;,
$$
.. math::
:nowrap:
$$
s_\nu = \frac{(1+\phi)}{2\tau_n} +\frac{(1-\phi)}{2\tau_w} \;,
$$
where
.. math::
:nowrap:
$$
\nu_w = \frac{1}{3}\Big(\tau_w - \frac{1}{2} \Big) \;, \quad
\nu_n = \frac{1}{3}\Big(\tau_n - \frac{1}{2} \Big) \;.
$$
These values are then used to model the momentum transport.
A D3Q19 LBE is constructed to describe the momentum transport
.. math::
:nowrap:
$$
f_q(\bm{x}_i + \bm{\xi}_q \delta t,t + \delta t) - f_q(\bm{x}_i,t) =
\sum^{Q-1}_{k=0} M^{-1}_{qk} S_{kk} (m_k^{eq}-m_k) + \sum^{Q-1}_{k=0} M^{-1}_{qk} (1-\frac{S_{kk}}{2}) \hat{F}_q\;,
$$
The force is imposed based on the construction developed by Guo et al
.. math::
:nowrap:
$$
F_i = \rho_0 \omega_i \left[\frac{\bm{e}_i \cdot \bm{a}}{c_s^2} +
\frac{\bm{u} \bm{a}:(\bm{e}_i \bm{e}_i -c_s^2 \mathcal{I})}{\epsilon c_s^4} \right] ,
$$
The acceleration includes contributions due to the external driving force :math:`\bm{g}`
as well as a drag force due to the permeability :math:`K` and flow rate :math:`\bm{u}` with the
porosity :math:`\epsilon` and viscosity :math:`\nu` determining the net forces acting within
a grey voxel
.. math::
:nowrap:
$$
\bm{a} = - \frac{\epsilon \nu}{K} \bm{u} + \bm{F}_{cp}/\rho_0 + \epsilon \bm{g},
$$
The flow velocity is defined as
.. math::
:nowrap:
$$
\rho_0 \bm{u} = \sum_i \bm{e}_i f_i + \frac{\delta t}{2} \rho_0 \bm{a}.
$$
Combining the previous expressions,
.. math::
:nowrap:
$$
\bm{u} = \frac{\frac{1}{\rho_0}\sum_i \bm{e}_i f_i + \frac{\delta t}{2}\epsilon \bm{g} +
\frac{\delta t}{2} \frac{\bm{F}_{cp}}{\rho_0}}{1+ \frac{\delta t}{2} \frac{\epsilon \nu}{K}}
$$
Where :math:`\bm{F}` is an external body force and :math:`c_s^2 = 1/3` is the speed of sound for the LB model.
The moments are linearly indepdendent:
.. math::
:nowrap:
$$
m_k = \sum_{q=0}^{18} M_{qk} f_q\;.
$$
The relaxation parameters are determined from the relaxation time:
.. math::
:nowrap:
$$
\lambda_1 = \lambda_2= \lambda_9 = \lambda_{10}= \lambda_{11}= \lambda_{12}= \lambda_{13}= \lambda_{14}= \lambda_{15} = s_\nu \;,
$$
.. math::
:nowrap:
$$
\lambda_{4}= \lambda_{6}= \lambda_{8} = \lambda_{16} = \lambda_{17} = \lambda_{18}= \frac{8(2-s_\nu)}{8-s_\nu} \;,
$$
The non-zero equilibrium moments are defined as
.. math::
:nowrap:
$$
m_1^{eq} = (j_x^2+j_y^2+j_z^2) - \alpha |\textbf{C}|, \\
$$
.. math::
:nowrap:
$$
m_9^{eq} = (2j_x^2-j_y^2-j_z^2)+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{11}^{eq} = (j_y^2-j_z^2) + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{13}^{eq} = j_x j_y + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\
$$
.. math::
:nowrap:
$$
m_{14}^{eq} = j_y j_z + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\
$$
.. math::
:nowrap:
$$
m_{15}^{eq} = j_x j_z + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;,
$$
where the color gradient is determined from the phase indicator field
.. math::
:nowrap:
$$
\textbf{C}=\nabla \phi\;.
$$
and the unit normal vector is
.. math::
:nowrap:
$$
\bm{n} = \frac{\textbf{C}}{|\textbf{C}|}\;.
$$
****************************
Boundary Conditions
****************************
Due to the nature of the contribution of the porosity to the pressure term in the Chapman-Enskog
expansion, periodic boundary conditions are recommended for ``lbpm_greyscaleColor_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
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

2
docs/source/userGuide/models/index.rst
View File

@ -18,6 +18,8 @@ Currently supported lattice Boltzmann models
greyscale/*
greyscaleColor/*
freeEnergy/*

35
docs/source/userGuide/models/mrt/mrt.rst
View File

@ -22,7 +22,7 @@ For additional details please refer to your local system documentation.
Model parameters
***************************
The essential model parameters for the color model are
The essential model parameters for the single-phase MRT model are
- ``tau`` -- control the fluid viscosity -- :math:`0.7 < \tau < 1.5`
@ -165,13 +165,6 @@ The relaxation parameters are determined based on the relaxation time :math:`\ta
$$
****************************
Example Input File
****************************
****************************
Boundary Conditions
@ -207,3 +200,29 @@ the inlet or outlet, the ``Domain`` section of the database may specify the foll
- ``InletLayerPhase = 2`` -- establish a reservoir of component B at the inlet
- ``OutletLayerPhase = 1`` -- establish a reservoir of component A at the outlet
****************************
Example Input File
****************************
.. code-block:: c
MRT {
tau = 1.0
F = 0.0, 0.0, 1.0e-5
timestepMax = 2000
tolerance = 0.01
}
Domain {
Filename = "Bentheimer_LB_sim_intermediate_oil_wet_Sw_0p37.raw"
ReadType = "8bit" // data type
N = 900, 900, 1600 // size of original image
nproc = 2, 2, 2 // process grid
n = 200, 200, 200 // sub-domain size
offset = 300, 300, 300 // offset to read sub-domain
voxel_length = 1.66 // voxel length (in microns)
ReadValues = 0, 1, 2 // labels within the original image
WriteValues = 0, 1, 2 // associated labels to be used by LBPM
InletLayers = 0, 0, 10 // specify 10 layers along the z-inlet
BC = 0 // boundary condition type (0 for periodic)
}

46
docs/source/userGuide/visualization/visit.rst
View File

@ -2,4 +2,48 @@
Visualizing simulation data with visit
======================================
placeholder for visit
Control over visualization options is set within the ``Visualization`` section of the input file
.. code:: c
Visualization {
write_silo = true // write SILO databases with assigned variables
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
save_phase_field = true // save phase field within SILO database
save_pressure = false // save pressure field within SILO database
save_velocity = false // save velocity field within SILO database
}
LBPM provides two main options for visualization. The first is the writing of 8-bit raw binary files, which are labeled based on the timestep. For example, if ``visualization_interval = 10000`` (specified within the Analysis section of the input file) the first 8-bit binary file will be written when ``timestep = 1000`` and will be named ``id_t1000.raw``. Additional files will be written subsequently at the specified interval. Similarly, higher fidelity visualization files are written using the SILO format, which are stored within the directories ``vis1000/``. The summary file ``LBM.visit`` enumerates these files so that they can be loaded directly into VisIt or other visualization software. By default, only the phase field will be saved. Visualization for other variables, such as the pressure and velocity fields, can be enabled by setting the associated flags to ``true``.
The VisIt software is able to natively read the SILO format. To import the data fields written by LBPM, open the VisIt GUI and select ``File > Open file`` from the top menu. Then select the LBM.visit file that you would like to read
.. figure:: ../../_static/images/lbpm-visit-workflow-i.png
:width: 600
:alt: VisIt GUI
Opening data in the VisIt GUI.
Once the file has been opened, any database fields within the file will be visible from within the program. Here we construct an 3D countour the phase field to visualize the boundary of the oil. The menu for the ``Contour`` object will show that the data for phase is available. The pressure and velocity fields will only be visible if they have been explicitly enabled within the simulation options (see ``Visualization`` details above)
.. figure:: ../../_static/images/lbpm-visit-workflow-ii.png
:width: 600
:alt: VisIt GUI
Selecting isosurface contour to represent a surface.
Once the contour has been created, double click the object icon to open the countour plot attributes window. In the field for ``Select by`` choose ``values`` and specify ``0`` as the contour to draw. Note that it is possible to change the color and opacity of the contour, which is useful if you want to draw multiple contours on the same plot
.. figure:: ../../_static/images/lbpm-visit-workflow-iii.png
:width: 600
:alt: VisIt GUI
Drawing an isosurface.
Once the attributes have been selected, click the Draw button to render the contour. Depending on the machine where you are rendering and the size of the image, it may take several minutes to render the window
.. figure:: ../../_static/images/lbpm-visit-workflow-iv.png
:width: 600
:alt: VisIt GUI
Rendering an isosurface.

1
doxygen/DoxygenMainpage.h
View File

@ -3,6 +3,7 @@
* C/C++ routines
*
* - \ref ScaLBL.h "Scalable Lattice Boltzmann Library (ScaLBL)"
* - \ref Domain.h "Domain structure"
* - \ref models "Lattice Boltzmann models"
* - \ref ScaLBL_ColorModel "Color model"
* - \ref analysis "Analysis routines"

6
example/DiscPack/discs.csv Normal file
View File

@ -0,0 +1,6 @@
cx cy r
0.25 0.2 0.12
0.45 0.52 0.15
0.7 0.8 0.12
0.8 0.25 0.13
0.2 0.85 0.1
1 cx cy r
2 0.25 0.2 0.12
3 0.45 0.52 0.15
4 0.7 0.8 0.12
5 0.8 0.25 0.13
6 0.2 0.85 0.1

42
hip/Ion.cu
View File

@ -98,13 +98,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *D
}
}
__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
__global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
double uEPx,uEPy,uEPz;//electrochemical induced velocity
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
int nr1,nr2,nr3,nr4,nr5,nr6;
@ -147,6 +148,20 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
nr6 = neighborList[n+5*Np];
f6 = dist[nr6];
// compute diffusive flux
flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
//dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
@ -178,13 +193,14 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
__global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
double uEPx,uEPy,uEPz;//electrochemical induced velocity
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
int S = Np/NBLOCKS/NTHREADS + 1;
@ -213,6 +229,20 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *V
f5 = dist[6*Np+n];
f6 = dist[5*Np+n];
// compute diffusive flux
flux_diffusive_x = (1.0-0.5*rlx)*((f1-f2)-ux*Ci);
flux_diffusive_y = (1.0-0.5*rlx)*((f3-f4)-uy*Ci);
flux_diffusive_z = (1.0-0.5*rlx)*((f5-f6)-uz*Ci);
FluxDiffusive[n+0*Np] = flux_diffusive_x;
FluxDiffusive[n+1*Np] = flux_diffusive_y;
FluxDiffusive[n+2*Np] = flux_diffusive_z;
FluxAdvective[n+0*Np] = ux*Ci;
FluxAdvective[n+1*Np] = uy*Ci;
FluxAdvective[n+2*Np] = uz*Ci;
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
//dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
@ -331,10 +361,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAodd_Ion<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Den,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
dvc_ScaLBL_D3Q7_AAodd_Ion<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
@ -343,10 +373,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAeven_Ion<<<NBLOCKS,NTHREADS >>>(dist,Den,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
dvc_ScaLBL_D3Q7_AAeven_Ion<<<NBLOCKS,NTHREADS >>>(dist,Den,FluxDiffusive,FluxAdvective,FluxElectrical,Velocity,ElectricField,Di,zi,rlx,Vt,start,finish,Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){

1
models/ColorModel.h
View File

@ -111,7 +111,6 @@ public:
std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
//std::shared_ptr<TwoPhase> Averages;
std::shared_ptr<SubPhase> Averages;
// input database

785
models/GreyscaleColorModel.cpp
View File

@ -356,8 +356,8 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
AFFINITY=AffinityList[idx];
Sn = SnList[idx];
Sw = SwList[idx];
Kn = SnList[idx];
Kw = SwList[idx];
Kn = KnList[idx];
Kw = KwList[idx];
idx = NLABELS;
}
}
@ -523,71 +523,6 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
delete [] Permeability;
}
//void ScaLBL_GreyscaleColorModel::AssignGreyscalePotential()
//{
// double *psi;//greyscale potential
// psi = new double[N];
//
// size_t NLABELS=0;
// signed char VALUE=0;
// double AFFINITY=0.f;
//
// auto LabelList = greyscaleColor_db->getVector<int>( "ComponentLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "ComponentAffinity" );
// NLABELS=LabelList.size();
//
// //first, copy over normal phase field
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// idx = NLABELS;
// }
// }
// // fluid labels are reserved
// if (VALUE == 1) AFFINITY=1.0;
// else if (VALUE == 2) AFFINITY=-1.0;
// psi[n] = AFFINITY;
// }
// }
// }
//
// //second, scale the phase field for grey nodes
// double Cap_Penalty=1.f;
// auto GreyLabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto PermeabilityList = greyscaleColor_db->getVector<double>( "PermeabilityList" );
// NLABELS=GreyLabelList.size();
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// Cap_Penalty=1.f;
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// if (VALUE == GreyLabelList[idx]){
// Cap_Penalty=alpha*W/sqrt(PermeabilityList[idx]/Dm->voxel_length/Dm->voxel_length);
// idx = NLABELS;
// }
// }
// //update greyscale potential
// psi[n] = psi[n]*Cap_Penalty;
// }
// }
// }
//
// ScaLBL_CopyToDevice(Psi, psi, N*sizeof(double));
// ScaLBL_Comm->Barrier();
// delete [] psi;
//}
void ScaLBL_GreyscaleColorModel::Create(){
/*
* This function creates the variables needed to run a LBM
@ -635,7 +570,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
//ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz);//greyscale potential
ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
//ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &MobilityRatio, sizeof(double)*Np);
//ScaLBL_AllocateDeviceMemory((void **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz);
//ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreySolidW, sizeof(double)*Np);
@ -686,8 +621,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
AssignComponentLabels();//do open/black/grey nodes initialization
AssignGreySolidLabels();
AssignGreyPoroPermLabels();
//AssignGreyscalePotential();
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity);
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
}
@ -787,9 +721,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
int nprocs=nprocx*nprocy*nprocz;
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
int IMAGE_INDEX = 0;
int IMAGE_COUNT = 0;
std::vector<std::string> ImageList;
bool SET_CAPILLARY_NUMBER = false;
bool RESCALE_FORCE = false;
bool MORPH_ADAPT = false;
@ -845,16 +776,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
/* defaults for simulation protocols */
auto protocol = greyscaleColor_db->getWithDefault<std::string>( "protocol", "none" );
if (protocol == "image sequence"){
// Get the list of images
USE_DIRECT = true;
ImageList = greyscaleColor_db->getVector<std::string>( "image_sequence");
IMAGE_INDEX = greyscaleColor_db->getWithDefault<int>( "image_index", 0 );
IMAGE_COUNT = ImageList.size();
morph_interval = 10000;
USE_MORPH = true;
}
else if (protocol == "seed water"){
if (protocol == "seed water"){
morph_delta = -0.05;
seed_water = 0.01;
USE_SEED = true;
@ -908,15 +830,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
if (rank==0){
printf("********************************************************\n");
if (protocol == "image sequence"){
printf(" using protocol = image sequence \n");
printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
printf(" tolerance = %f \n",tolerance);
std::string first_image = ImageList[IMAGE_INDEX];
printf(" first image in sequence: %s ***\n", first_image.c_str());
}
else if (protocol == "seed water"){
if (protocol == "seed water"){
printf(" using protocol = seed water \n");
printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
@ -963,18 +877,9 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
// Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,MobilityRatio,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//Model-1&4
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
////Model-2&3
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
@ -992,18 +897,9 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,MobilityRatio,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
//Model-1&4
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
////Model-2&3
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier();
// *************EVEN TIMESTEP*************
@ -1025,18 +921,9 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,MobilityRatio,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//Model-1&4
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
////Model-2&3
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
@ -1054,18 +941,9 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,MobilityRatio,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
//Model-1&4
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
////Model-2&3
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier();
//************************************************************************
PROFILE_STOP("Update");
@ -1121,11 +999,13 @@ void ScaLBL_GreyscaleColorModel::Run(){
if (timestep%analysis_interval == 0){
ScaLBL_Comm->RegularLayout(Map,Pressure,Averages->Pressure);
ScaLBL_Comm->RegularLayout(Map,MobilityRatio,Averages->MobilityRatio);
ScaLBL_Comm->RegularLayout(Map,&Den[0],Averages->Rho_n);
ScaLBL_Comm->RegularLayout(Map,&Den[Np],Averages->Rho_w);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Averages->Vel_x);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Averages->Vel_y);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Averages->Vel_z);
Averages->Basic();
}
@ -1212,43 +1092,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
double pA = Averages->Oil.p;
double pB = Averages->Water.p;
double pAB = (pA-pB)/(h*6.0*alpha);
// -------- The following quantities may not make sense for greyscale simulation -----------//
// double pAc = Averages->gnc.p;
// double pBc = Averages->gwc.p;
// double pAB_connected = (pAc-pBc)/(h*6.0*alpha);
// // connected contribution
// double Vol_nc = Averages->gnc.V/Dm->Volume;
// double Vol_wc = Averages->gwc.V/Dm->Volume;
// double Vol_nd = Averages->gnd.V/Dm->Volume;
// double Vol_wd = Averages->gwd.V/Dm->Volume;
// double Mass_n = Averages->gnc.M + Averages->gnd.M;
// double Mass_w = Averages->gwc.M + Averages->gwd.M;
// double vAc_x = Averages->gnc.Px/Mass_n;
// double vAc_y = Averages->gnc.Py/Mass_n;
// double vAc_z = Averages->gnc.Pz/Mass_n;
// double vBc_x = Averages->gwc.Px/Mass_w;
// double vBc_y = Averages->gwc.Py/Mass_w;
// double vBc_z = Averages->gwc.Pz/Mass_w;
// // disconnected contribution
// double vAd_x = Averages->gnd.Px/Mass_n;
// double vAd_y = Averages->gnd.Py/Mass_n;
// double vAd_z = Averages->gnd.Pz/Mass_n;
// double vBd_x = Averages->gwd.Px/Mass_w;
// double vBd_y = Averages->gwd.Py/Mass_w;
// double vBd_z = Averages->gwd.Pz/Mass_w;
//
// double flow_rate_A_connected = Vol_nc*(vAc_x*dir_x + vAc_y*dir_y + vAc_z*dir_z);
// double flow_rate_B_connected = Vol_wc*(vBc_x*dir_x + vBc_y*dir_y + vBc_z*dir_z);
// double flow_rate_A_disconnected = (Vol_nd)*(vAd_x*dir_x + vAd_y*dir_y + vAd_z*dir_z);
// double flow_rate_B_disconnected = (Vol_wd)*(vBd_x*dir_x + vBd_y*dir_y + vBd_z*dir_z);
//
// double kAeff_connected = h*h*muA*flow_rate_A_connected/(force_mag);
// double kBeff_connected = h*h*muB*flow_rate_B_connected/(force_mag);
//
// double kAeff_disconnected = h*h*muA*flow_rate_A_disconnected/(force_mag);
// double kBeff_disconnected = h*h*muB*flow_rate_B_disconnected/(force_mag);
// //---------------------------------------------------------------------------------------//
double kAeff = h*h*muA*(flow_rate_A)/(force_mag);
double kBeff = h*h*muB*(flow_rate_B)/(force_mag);
@ -1300,22 +1143,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
if (MORPH_ADAPT ){
CURRENT_MORPH_TIMESTEPS += analysis_interval;
if (USE_DIRECT){
// Use image sequence
IMAGE_INDEX++;
MORPH_ADAPT = false;
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);
greyscaleColor_db->putScalar<int>("image_index",IMAGE_INDEX);
ImageInit(next_image);
}
else{
if (rank==0) printf("Finished simulating image sequence \n");
timestep = timestepMax;
}
}
else if (USE_SEED){
if (USE_SEED){
delta_volume = volA*Dm->Volume - initial_volume;
CURRENT_MORPH_TIMESTEPS += analysis_interval;
double massChange = SeedPhaseField(seed_water);
@ -1366,50 +1194,6 @@ void ScaLBL_GreyscaleColorModel::Run(){
// ************************************************************************
}
void ScaLBL_GreyscaleColorModel::ImageInit(std::string Filename){
if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str());
Mask->Decomp(Filename);
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; // save what was read
AssignComponentLabels();
AssignGreySolidLabels();
AssignGreyPoroPermLabels();
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity);
ScaLBL_Comm->RegularLayout(Map,Porosity_dvc,Averages->Porosity);
//NOTE in greyscale simulations, water may have multiple labels (e.g. 2, 21, 22, etc)
//so the saturaiton calculation is not that straightforward
// double Count = 0.0;
// double PoreCount = 0.0;
// for (int k=1; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
// if (id[Nx*Ny*k+Nx*j+i] == 2){
// PoreCount++;
// Count++;
// }
// else if (id[Nx*Ny*k+Nx*j+i] == 1){
// PoreCount++;
// }
// }
// }
// }
// Count=Dm->Comm.sumReduce( Count);
// PoreCount=Dm->Comm.sumReduce( PoreCount);
// if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
ScaLBL_D3Q19_Init(fq, Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->Barrier();
//ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double));
//double saturation = Count/PoreCount;
//return saturation;
}
double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil){
srand(time(NULL));
double mass_loss =0.f;
@ -1713,546 +1497,3 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
*/
}
//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-1
//{
// // ONLY initialize grey nodes
// // Key input parameters:
// // 1. GreySolidLabels
// // labels for grey nodes
// // 2. GreySolidAffinity
// // affinity ranges [-1,1]
// // oil-wet > 0
// // water-wet < 0
// // neutral = 0
// double *SolidPotential_host = new double [Nx*Ny*Nz];
// double *GreySolidGrad_host = new double [3*Np];
//
// size_t NLABELS=0;
// signed char VALUE=0;
// double AFFINITY=0.f;
//
// auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
//
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
// }
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// idx = NLABELS;
// //Mask->id[n] = 0; // set mask to zero since this is an immobile component
// }
// }
// SolidPotential_host[n] = AFFINITY;
// }
// }
// }
//
// // Calculate grey-solid color-gradient
// double *Dst;
// Dst = new double [3*3*3];
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
// int index = kk*9+jj*3+ii;
// Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
// }
// }
// }
// double w_face = 1.f;
// double w_edge = 0.5;
// double w_corner = 0.f;
// //local
// Dst[13] = 0.f;
// //faces
// Dst[4] = w_face;
// Dst[10] = w_face;
// Dst[12] = w_face;
// Dst[14] = w_face;
// Dst[16] = w_face;
// Dst[22] = w_face;
// // corners
// Dst[0] = w_corner;
// Dst[2] = w_corner;
// Dst[6] = w_corner;
// Dst[8] = w_corner;
// Dst[18] = w_corner;
// Dst[20] = w_corner;
// Dst[24] = w_corner;
// Dst[26] = w_corner;
// // edges
// Dst[1] = w_edge;
// Dst[3] = w_edge;
// Dst[5] = w_edge;
// Dst[7] = w_edge;
// Dst[9] = w_edge;
// Dst[11] = w_edge;
// Dst[15] = w_edge;
// Dst[17] = w_edge;
// Dst[19] = w_edge;
// Dst[21] = w_edge;
// Dst[23] = w_edge;
// Dst[25] = w_edge;
//
// for (int k=1; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
// int idx=Map(i,j,k);
// if (!(idx < 0)){
// double phi_x = 0.f;
// double phi_y = 0.f;
// double phi_z = 0.f;
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
//
// int index = kk*9+jj*3+ii;
// double weight= Dst[index];
//
// int idi=i+ii-1;
// int idj=j+jj-1;
// int idk=k+kk-1;
//
// if (idi < 0) idi=0;
// if (idj < 0) idj=0;
// if (idk < 0) idk=0;
// if (!(idi < Nx)) idi=Nx-1;
// if (!(idj < Ny)) idj=Ny-1;
// if (!(idk < Nz)) idk=Nz-1;
//
// int nn = idk*Nx*Ny + idj*Nx + idi;
// double vec_x = double(ii-1);
// double vec_y = double(jj-1);
// double vec_z = double(kk-1);
// double GWNS=SolidPotential_host[nn];
// phi_x += GWNS*weight*vec_x;
// phi_y += GWNS*weight*vec_y;
// phi_z += GWNS*weight*vec_z;
// }
// }
// }
// GreySolidGrad_host[idx+0*Np] = phi_x;
// GreySolidGrad_host[idx+1*Np] = phi_y;
// GreySolidGrad_host[idx+2*Np] = phi_z;
// }
// }
// }
// }
//
// if (rank==0){
// printf("Number of Grey-solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// VALUE=LabelList[idx];
// AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
//
// ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
// ScaLBL_Comm->Barrier();
// delete [] SolidPotential_host;
// delete [] GreySolidGrad_host;
// delete [] Dst;
//}
////----------------------------------------------------------------------------------------------------------//
//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-2 & Model-3
//{
// // ONLY initialize grey nodes
// // Key input parameters:
// // 1. GreySolidLabels
// // labels for grey nodes
// // 2. GreySolidAffinity
// // affinity ranges [-1,1]
// // oil-wet > 0
// // water-wet < 0
// // neutral = 0
//
// double *GreySolidPhi_host = new double [Nx*Ny*Nz];
// //initialize grey solid 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;
// GreySolidPhi_host[n]=0.f;
// }
// }
// }
//
// auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
//
// size_t NLABELS=0;
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
// }
//
// double *Dst;
// Dst = new double [3*3*3];
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
// int index = kk*9+jj*3+ii;
// Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
// }
// }
// }
// double w_face = 1.f;
// double w_edge = 1.f;
// double w_corner = 0.f;
// //local
// Dst[13] = 0.f;
// //faces
// Dst[4] = w_face;
// Dst[10] = w_face;
// Dst[12] = w_face;
// Dst[14] = w_face;
// Dst[16] = w_face;
// Dst[22] = w_face;
// // corners
// Dst[0] = w_corner;
// Dst[2] = w_corner;
// Dst[6] = w_corner;
// Dst[8] = w_corner;
// Dst[18] = w_corner;
// Dst[20] = w_corner;
// Dst[24] = w_corner;
// Dst[26] = w_corner;
// // edges
// Dst[1] = w_edge;
// Dst[3] = w_edge;
// Dst[5] = w_edge;
// Dst[7] = w_edge;
// Dst[9] = w_edge;
// Dst[11] = w_edge;
// Dst[15] = w_edge;
// Dst[17] = w_edge;
// Dst[19] = w_edge;
// Dst[21] = w_edge;
// Dst[23] = w_edge;
// Dst[25] = w_edge;
//
// for (int k=1; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
//
// int n = k*Nx*Ny+j*Nx+i;
// signed char VALUE=Mask->id[n];
// double AFFINITY=0.f;
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// idx = NLABELS;
// //Mask->id[n] = 0; // set mask to zero since this is an immobile component
// }
// }
//
// if (VALUE>2){//i.e. a grey node
// double neighbor_counter = 0;
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
//
// int index = kk*9+jj*3+ii;
// double weight= Dst[index];
//
// int idi=i+ii-1;
// int idj=j+jj-1;
// int idk=k+kk-1;
//
// if (idi < 0) idi=0;
// if (idj < 0) idj=0;
// if (idk < 0) idk=0;
// if (!(idi < Nx)) idi=Nx-1;
// if (!(idj < Ny)) idj=Ny-1;
// if (!(idk < Nz)) idk=Nz-1;
//
// int nn = idk*Nx*Ny + idj*Nx + idi;
// //if (Mask->id[nn] != VALUE){//Model-2:i.e. open nodes, impermeable solid nodes or any other type of greynodes
// if (Mask->id[nn] <=0){//Model-3:i.e. only impermeable solid nodes or any other type of greynodes
// neighbor_counter +=weight;
// }
// }
// }
// }
// if (neighbor_counter>0){
// GreySolidPhi_host[n] = AFFINITY;
// }
// }
// }
// }
// }
//
// if (rank==0){
// printf("Number of grey-solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// signed char VALUE=LabelList[idx];
// double AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
// ScaLBL_CopyToDevice(GreySolidPhi, GreySolidPhi_host, Nx*Ny*Nz*sizeof(double));
// ScaLBL_Comm->Barrier();
//
// //debug
// //FILE *OUTFILE;
// //sprintf(LocalRankFilename,"GreySolidInit.%05i.raw",rank);
// //OUTFILE = fopen(LocalRankFilename,"wb");
// //fwrite(GreySolidPhi_host,8,N,OUTFILE);
// //fclose(OUTFILE);
//
// delete [] GreySolidPhi_host;
// delete [] Dst;
//}
//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
//{
// // ONLY initialize grey nodes
// // Key input parameters:
// // 1. GreySolidLabels
// // labels for grey nodes
// // 2. GreySolidAffinity
// // affinity ranges [-1,1]
// // oil-wet > 0
// // water-wet < 0
// // neutral = 0
// double *SolidPotential_host = new double [Nx*Ny*Nz];
// double *GreySolidGrad_host = new double [3*Np];
//
// size_t NLABELS=0;
// signed char VALUE=0;
// double AFFINITY=0.f;
//
// auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
//
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
// }
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// idx = NLABELS;
// //Mask->id[n] = 0; // set mask to zero since this is an immobile component
// }
// }
// SolidPotential_host[n] = AFFINITY;
// }
// }
// }
//
// // Calculate grey-solid color-gradient
// double *Dst;
// Dst = new double [3*3*3];
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
// int index = kk*9+jj*3+ii;
// Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
// }
// }
// }
// double w_face = 1.f;
// double w_edge = 0.5;
// double w_corner = 0.f;
// //local
// Dst[13] = 0.f;
// //faces
// Dst[4] = w_face;
// Dst[10] = w_face;
// Dst[12] = w_face;
// Dst[14] = w_face;
// Dst[16] = w_face;
// Dst[22] = w_face;
// // corners
// Dst[0] = w_corner;
// Dst[2] = w_corner;
// Dst[6] = w_corner;
// Dst[8] = w_corner;
// Dst[18] = w_corner;
// Dst[20] = w_corner;
// Dst[24] = w_corner;
// Dst[26] = w_corner;
// // edges
// Dst[1] = w_edge;
// Dst[3] = w_edge;
// Dst[5] = w_edge;
// Dst[7] = w_edge;
// Dst[9] = w_edge;
// Dst[11] = w_edge;
// Dst[15] = w_edge;
// Dst[17] = w_edge;
// Dst[19] = w_edge;
// Dst[21] = w_edge;
// Dst[23] = w_edge;
// Dst[25] = w_edge;
//
// for (int k=1; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
// int idx=Map(i,j,k);
// if (!(idx < 0)){
// double phi_x = 0.f;
// double phi_y = 0.f;
// double phi_z = 0.f;
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
//
// int index = kk*9+jj*3+ii;
// double weight= Dst[index];
//
// int idi=i+ii-1;
// int idj=j+jj-1;
// int idk=k+kk-1;
//
// if (idi < 0) idi=0;
// if (idj < 0) idj=0;
// if (idk < 0) idk=0;
// if (!(idi < Nx)) idi=Nx-1;
// if (!(idj < Ny)) idj=Ny-1;
// if (!(idk < Nz)) idk=Nz-1;
//
// int nn = idk*Nx*Ny + idj*Nx + idi;
// double vec_x = double(ii-1);
// double vec_y = double(jj-1);
// double vec_z = double(kk-1);
// double GWNS=SolidPotential_host[nn];
// phi_x += GWNS*weight*vec_x;
// phi_y += GWNS*weight*vec_y;
// phi_z += GWNS*weight*vec_z;
// }
// }
// }
// if (Averages->SDs(i,j,k)<2.0){
// GreySolidGrad_host[idx+0*Np] = phi_x;
// GreySolidGrad_host[idx+1*Np] = phi_y;
// GreySolidGrad_host[idx+2*Np] = phi_z;
// }
// else{
// GreySolidGrad_host[idx+0*Np] = 0.0;
// GreySolidGrad_host[idx+1*Np] = 0.0;
// GreySolidGrad_host[idx+2*Np] = 0.0;
// }
// }
// }
// }
// }
//
//
// if (rank==0){
// printf("Number of Grey-solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// VALUE=LabelList[idx];
// AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
//
// ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
// ScaLBL_Comm->Barrier();
// delete [] SolidPotential_host;
// delete [] GreySolidGrad_host;
// delete [] Dst;
//}
//--------- This is another old version of calculating greyscale-solid color-gradient modification-------//
// **not working effectively, to be deprecated
//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()
//{
// // ONLY initialize grey nodes
// // Key input parameters:
// // 1. GreySolidLabels
// // labels for grey nodes
// // 2. GreySolidAffinity
// // affinity ranges [-1,1]
// // oil-wet > 0
// // water-wet < 0
// // neutral = 0
//
// //double *SolidPotential_host = new double [Nx*Ny*Nz];
// double *GreySolidPhi_host = new double [Nx*Ny*Nz];
// signed char VALUE=0;
// double AFFINITY=0.f;
//
// auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
//
// size_t NLABELS=0;
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
// }
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// idx = NLABELS;
// //Mask->id[n] = 0; // set mask to zero since this is an immobile component
// }
// }
// GreySolidPhi_host[n] = AFFINITY;
// }
// }
// }
//
// if (rank==0){
// printf("Number of grey-solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// VALUE=LabelList[idx];
// AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, solid-affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
// ScaLBL_CopyToDevice(GreySolidPhi, GreySolidPhi_host, Nx*Ny*Nz*sizeof(double));
// ScaLBL_Comm->Barrier();
//
// //debug
// FILE *OUTFILE;
// sprintf(LocalRankFilename,"GreySolidInit.%05i.raw",rank);
// OUTFILE = fopen(LocalRankFilename,"wb");
// fwrite(GreySolidPhi_host,8,N,OUTFILE);
// fclose(OUTFILE);
//
// //delete [] SolidPotential_host;
// delete [] GreySolidPhi_host;
//}
//----------------------------------------------------------------------------------------------------------//

72
models/GreyscaleColorModel.h
View File

@ -15,20 +15,71 @@ Implementation of two-fluid greyscale color lattice boltzmann model
#include "ProfilerApp.h"
#include "threadpool/thread_pool.h"
/**
* \class ScaLBL_GreyscaleColorModel
*
* @details
* The ScaLBL_GreyscaleColorModel class extends the standard color model incorporate transport
* through sub-resolution "greyscale" regions.
* Momentum transport equations are described by a D3Q19 scheme
* Mass transport equations are described by D3Q7 scheme
*/
class ScaLBL_GreyscaleColorModel{
public:
/**
* \brief Constructor
* @param RANK processor rank
* @param NP number of processors
* @param COMM MPI communicator
*/
ScaLBL_GreyscaleColorModel(int RANK, int NP, const Utilities::MPI& COMM);
~ScaLBL_GreyscaleColorModel();
// functions in they should be run
/**
* \brief Read simulation parameters
* @param filename input database file that includes "Color" section
*/
void ReadParams(string filename);
/**
* \brief Read simulation parameters
* @param db0 input database that includes "Color" section
*/
void ReadParams(std::shared_ptr<Database> db0);
/**
* \brief Create domain data structures
*/
void SetDomain();
/**
* \brief Read image data
*/
void ReadInput();
/**
* \brief Create color model data structures
*/
void Create();
/**
* \brief Initialize the simulation
*/
void Initialize();
/**
* \brief Run the simulation
*/
void Run();
/**
* \brief Debugging function to dump simulation state to disk
*/
void WriteDebug();
void WriteVisFiles();
bool Restart,pBC;
bool REVERSE_FLOW_DIRECTION;
@ -72,7 +123,7 @@ public:
double *GreySw;
double *GreyKn;
double *GreyKw;
//double *ColorGrad;
double *MobilityRatio;
double *Velocity;
double *Pressure;
double *Porosity_dvc;
@ -91,14 +142,23 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
/**
* \brief Assign wetting affinity values
*/
void AssignComponentLabels();
/**
* \brief Assign wetting affinity values in greyscale regions
*/
void AssignGreySolidLabels();
/**
* \brief Assign porosity and permeability in greyscale regions
*/
void AssignGreyPoroPermLabels();
//void AssignGreyscalePotential();
void ImageInit(std::string filename);
double MorphInit(const double beta, const double morph_delta);
/**
* \brief Seed phase field
*/
double SeedPhaseField(const double seed_water_in_oil);
double MorphOpenConnected(double target_volume_change);
void WriteVisFiles();
};

68
models/GreyscaleModel.cpp
View File

@ -617,29 +617,6 @@ void ScaLBL_GreyscaleModel::Run(){
if (BoundaryCondition > 0 && Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin = 1 + Dm->inlet_layers_z;//"1" indicates the halo layer
if (BoundaryCondition > 0 && Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax = Nz-1 - Dm->outlet_layers_z;
// px_loc = py_loc = pz_loc = 0.f;
// mass_loc = 0.f;
// for (int k=kmin; k<kmax; k++){
// for (int j=jmin; j<Ny-1; j++){
// for (int i=imin; i<Nx-1; i++){
// if (SignDist(i,j,k) > 0){
// px_loc += Velocity_x(i,j,k)*Den*PorosityMap(i,j,k);
// py_loc += Velocity_y(i,j,k)*Den*PorosityMap(i,j,k);
// pz_loc += Velocity_z(i,j,k)*Den*PorosityMap(i,j,k);
// mass_loc += Den*PorosityMap(i,j,k);
// }
// }
// }
// }
// MPI_Allreduce(&px_loc, &px, 1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
// MPI_Allreduce(&py_loc, &py, 1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
// MPI_Allreduce(&pz_loc, &pz, 1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
// MPI_Allreduce(&mass_loc,&mass_glb,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
//
// vax = px/mass_glb;
// vay = py/mass_glb;
// vaz = pz/mass_glb;
vax_loc = vay_loc = vaz_loc = 0.f;
for (int k=kmin; k<kmax; k++){
for (int j=jmin; j<Ny-1; j++){
@ -766,51 +743,6 @@ void ScaLBL_GreyscaleModel::Run(){
void ScaLBL_GreyscaleModel::VelocityField(){
/* Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); comm.barrier();
ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],3*SIZE);
memcpy(Morphology.SDn.data(), Distance.data(), Nx*Ny*Nz*sizeof(double));
Morphology.Initialize();
Morphology.UpdateMeshValues();
Morphology.ComputeLocal();
Morphology.Reduce();
double count_loc=0;
double count;
double vax,vay,vaz;
double vax_loc,vay_loc,vaz_loc;
vax_loc = vay_loc = vaz_loc = 0.f;
for (int n=0; n<ScaLBL_Comm->LastExterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
count_loc+=1.0;
}
for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
count_loc+=1.0;
}
MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
vax /= count;
vay /= count;
vaz /= count;
double mu = (tau-0.5)/3.f;
if (rank==0) printf("Fx Fy Fz mu Vs As Js Xs vx vy vz\n");
if (rank==0) printf("%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",Fx, Fy, Fz, mu,
Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz);
*/
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);

193
models/IonModel.cpp
View File

@ -692,6 +692,9 @@ void ScaLBL_IonModel::Create(){
ScaLBL_AllocateDeviceMemory((void **) &fq, number_ion_species*7*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Ci, number_ion_species*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &ChargeDensity, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &FluxDiffusive, number_ion_species*3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &FluxAdvective, number_ion_species*3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &FluxElectrical, number_ion_species*3*sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("LB Ion Solver: Setting up device map and neighbor list \n");
@ -877,9 +880,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//LB-Ion collison
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
if (BoundaryConditionSolid==1){
@ -933,9 +936,9 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//LB-Ion collison
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],&FluxDiffusive[3*ic*Np],&FluxAdvective[3*ic*Np],&FluxElectrical[3*ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
if (BoundaryConditionSolid==1){
@ -973,7 +976,54 @@ void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const s
ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration);
ScaLBL_Comm->Barrier(); comm.barrier();
IonConcentration_LB_to_Phys(IonConcentration);
}
void ScaLBL_IonModel::getIonFluxDiffusive(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic){
//This function wirte out the data in a normal layout (by aggregating all decomposed domains)
ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+0*Np],IonFlux_x);
IonFlux_LB_to_Phys(IonFlux_x,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+1*Np],IonFlux_y);
IonFlux_LB_to_Phys(IonFlux_y,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+2*Np],IonFlux_z);
IonFlux_LB_to_Phys(IonFlux_z,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
}
void ScaLBL_IonModel::getIonFluxAdvective(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic){
//This function wirte out the data in a normal layout (by aggregating all decomposed domains)
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+0*Np],IonFlux_x);
IonFlux_LB_to_Phys(IonFlux_x,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+1*Np],IonFlux_y);
IonFlux_LB_to_Phys(IonFlux_y,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+2*Np],IonFlux_z);
IonFlux_LB_to_Phys(IonFlux_z,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
}
void ScaLBL_IonModel::getIonFluxElectrical(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic){
//This function wirte out the data in a normal layout (by aggregating all decomposed domains)
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+0*Np],IonFlux_x);
IonFlux_LB_to_Phys(IonFlux_x,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+1*Np],IonFlux_y);
IonFlux_LB_to_Phys(IonFlux_y,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+2*Np],IonFlux_z);
IonFlux_LB_to_Phys(IonFlux_z,ic);
ScaLBL_Comm->Barrier(); comm.barrier();
}
void ScaLBL_IonModel::getIonConcentration_debug(int timestep){
@ -992,13 +1042,146 @@ void ScaLBL_IonModel::getIonConcentration_debug(int timestep){
}
}
void ScaLBL_IonModel::getIonFluxDiffusive_debug(int timestep){
//This function write out decomposed data
DoubleArray PhaseField(Nx,Ny,Nz);
for (size_t ic=0; ic<number_ion_species; ic++){
//x-component
ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+0*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_X;
sprintf(LocalRankFilename,"IonFluxDiffusive_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_X = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_X);
fclose(OUTFILE_X);
//y-component
ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+1*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Y;
sprintf(LocalRankFilename,"IonFluxDiffusive_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Y = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_Y);
fclose(OUTFILE_Y);
//z-component
ScaLBL_Comm->RegularLayout(Map,&FluxDiffusive[ic*3*Np+2*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Z;
sprintf(LocalRankFilename,"IonFluxDiffusive_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Z = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_Z);
fclose(OUTFILE_Z);
}
}
void ScaLBL_IonModel::getIonFluxAdvective_debug(int timestep){
//This function write out decomposed data
DoubleArray PhaseField(Nx,Ny,Nz);
for (size_t ic=0; ic<number_ion_species; ic++){
//x-component
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+0*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_X;
sprintf(LocalRankFilename,"IonFluxAdvective_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_X = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_X);
fclose(OUTFILE_X);
//y-component
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+1*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Y;
sprintf(LocalRankFilename,"IonFluxAdvective_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Y = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_Y);
fclose(OUTFILE_Y);
//z-component
ScaLBL_Comm->RegularLayout(Map,&FluxAdvective[ic*3*Np+2*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Z;
sprintf(LocalRankFilename,"IonFluxAdvective_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Z = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_Z);
fclose(OUTFILE_Z);
}
}
void ScaLBL_IonModel::getIonFluxElectrical_debug(int timestep){
//This function write out decomposed data
DoubleArray PhaseField(Nx,Ny,Nz);
for (size_t ic=0; ic<number_ion_species; ic++){
//x-component
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+0*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_X;
sprintf(LocalRankFilename,"IonFluxElectrical_X_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_X = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_X);
fclose(OUTFILE_X);
//y-component
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+1*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Y;
sprintf(LocalRankFilename,"IonFluxElectrical_Y_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Y = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_Y);
fclose(OUTFILE_Y);
//z-component
ScaLBL_Comm->RegularLayout(Map,&FluxElectrical[ic*3*Np+2*Np],PhaseField);
ScaLBL_Comm->Barrier(); comm.barrier();
IonFlux_LB_to_Phys(PhaseField,ic);
FILE *OUTFILE_Z;
sprintf(LocalRankFilename,"IonFluxElectrical_Z_%02zu_Time_%i.%05i.raw",ic+1,timestep,rank);
OUTFILE_Z = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE_Z);
fclose(OUTFILE_Z);
}
}
void ScaLBL_IonModel::IonConcentration_LB_to_Phys(DoubleArray &Den_reg){
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int idx=Map(i,j,k);
if (!(idx < 0)){
Den_reg(i,j,k) = Den_reg(i,j,k)/(h*h*h*1.0e-18);
Den_reg(i,j,k) = Den_reg(i,j,k)/(h*h*h*1.0e-18);//this converts the unit to [mol/m^3]
}
}
}
}
}
void ScaLBL_IonModel::IonFlux_LB_to_Phys(DoubleArray &Den_reg, const size_t ic){
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int idx=Map(i,j,k);
if (!(idx < 0)){
Den_reg(i,j,k) = Den_reg(i,j,k)/(h*h*1.0e-12*time_conv[ic]);//this converts the unit to [mol/m^2/s]
}
}
}

10
models/IonModel.h
View File

@ -36,6 +36,12 @@ public:
void Run(double *Velocity, double *ElectricField);
void getIonConcentration(DoubleArray &IonConcentration, const size_t ic);
void getIonConcentration_debug(int timestep);
void getIonFluxDiffusive(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic);
void getIonFluxAdvective(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic);
void getIonFluxElectrical(DoubleArray &IonFlux_x,DoubleArray &IonFlux_y,DoubleArray &IonFlux_z,const size_t ic);
void getIonFluxDiffusive_debug(int timestep);
void getIonFluxAdvective_debug(int timestep);
void getIonFluxElectrical_debug(int timestep);
void DummyFluidVelocity();
void DummyElectricField();
double CalIonDenConvergence(vector<double> &ci_avg_previous);
@ -83,6 +89,9 @@ public:
double *IonSolid;
double *FluidVelocityDummy;
double *ElectricFieldDummy;
double *FluxDiffusive;
double *FluxAdvective;
double *FluxElectrical;
private:
Utilities::MPI comm;
@ -98,5 +107,6 @@ private:
void AssignSolidBoundary(double *ion_solid);
void AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion,int ic);
void IonConcentration_LB_to_Phys(DoubleArray &Den_reg);
void IonFlux_LB_to_Phys(DoubleArray &Den_reg, const size_t ic);
};
#endif

11
tests/TestIonModel.cpp
View File

@ -21,13 +21,11 @@ int main(int argc, char **argv)
{
// Initialize MPI and error handlers
Utilities::startup( argc, argv );
Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank();
int nprocs = comm.getSize();
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
// Initialize MPI
Utilities::startup( argc, argv );
Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank();
int nprocs = comm.getSize();
if (rank == 0){
printf("**************************************\n");
@ -77,6 +75,9 @@ int main(int argc, char **argv)
}
}
IonModel.getIonConcentration_debug(timestep);
IonModel.getIonFluxDiffusive_debug(timestep);
IonModel.getIonFluxAdvective_debug(timestep);
IonModel.getIonFluxElectrical_debug(timestep);
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
if (rank==0) printf("*************************************************************\n");

2
tests/lbpm_greyscaleColor_simulator.cpp
View File

@ -53,7 +53,7 @@ int main(int argc, char **argv)
GreyscaleColor.Create(); // creating the model will create data structure to match the pore structure and allocate variables
GreyscaleColor.Initialize(); // initializing the model will set initial conditions for variables
GreyscaleColor.Run();
GreyscaleColor.WriteDebug();
GreyscaleColor.WriteVisFiles();
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_greyscaleColor_simulator",1);