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

merge doc update

Browse Source
This commit is contained in:
James McClure 2021-10-21 20:26:21 -04:00
commit 3f746e8ed2
50 changed files with 1808 additions and 2726 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);
}
}

1
analysis/FlowAdaptor.h
View File

@ -15,7 +15,6 @@
/**
* \class FlowAdaptor
*
* @brief
* The FlowAdaptor class operates on a lattice Boltzmann model to alter the flow conditions
*

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

@ -32,87 +32,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;
@ -140,6 +151,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; }
@ -181,22 +195,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];
@ -214,6 +284,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>

26
common/ScaLBL.h
View File

@ -15,7 +15,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
/** @file ScaLBL.h */
/* \detail Header file for Scalable Lattice Boltzmann Library
/* \details Header file for Scalable Lattice Boltzmann Library
* Separate implementations for GPU and CPU must both follow the conventions defined in this header
* This libarry contains the essential components of the LBM
* - streaming implementations
@ -48,7 +48,7 @@ extern "C" void ScaLBL_FreeDeviceMemory(void* pointer);
/**
* \brief Copy memory from host to device
* \detail Device memory should be close to simulation (based on NUMA cost)
* \details Device memory should be close to simulation (based on NUMA cost)
* Host memory may be a shared memory region (with possibly higher NUMA cost for simulation)
* Analysis routine should minimize NUMA for host memory (based on process placement)
* @param dest memory location to copy to
@ -60,7 +60,7 @@ extern "C" void ScaLBL_CopyToDevice(void* dest, const void* source, size_t size)
/**
* \brief Copy memory from device to host
* \detail Device memory should be close to simulation (based on NUMA cost)
* \details Device memory should be close to simulation (based on NUMA cost)
* Host memory may be a shared memory region (with possibly higher NUMA cost for simulation)
* Analysis routine should minimize NUMA for host memory (based on process placement)
* @param dest memory location to copy to
@ -78,7 +78,7 @@ extern "C" void ScaLBL_AllocateZeroCopy(void** address, size_t size);
/**
* \brief Copy memory from host to zero copy buffer
* \detail Device memory should be close to simulation (based on NUMA cost)
* \details Device memory should be close to simulation (based on NUMA cost)
* Host memory may be a shared memory region (with possibly higher NUMA cost for simulation)
* Analysis routine should minimize NUMA for host memory (based on process placement)
* @param dest memory location to copy to
@ -110,7 +110,7 @@ extern "C" void ScaLBL_D3Q19_Pack(int q, int *list, int start, int count, double
* @param list - list of distributions to communicate
* @param start - index to start parsing the list
* @param count - number of values to unppack
* @param revbuf - memory buffer where recieved values have been stored
* @param recvbuf - memory buffer where recieved values have been stored
* @param dist - memory buffer to hold the distributions
* @param N - size of the distributions (derived from Domain structure)
*/
@ -122,7 +122,7 @@ extern "C" void ScaLBL_D3Q19_Unpack(int q, int *list, int start, int count, doub
* @param list - list of distributions to communicate
* @param start - index to start parsing the list
* @param count - number of values to unppack
* @param revbuf - memory buffer where recieved values have been stored
* @param recvbuf - memory buffer where recieved values have been stored
* @param dist - memory buffer to hold the distributions
* @param N - size of the distributions (derived from Domain structure)
*/
@ -175,15 +175,15 @@ extern "C" void ScaLBL_D3Q19_Init(double *Dist, int Np);
/**
* \brief Compute momentum from D3Q19 distribution
* @param Dist - D3Q19 distributions
* @param dist - D3Q19 distributions
* @param vel - memory buffer to store the momentum that is computed
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np);
/**
* \brief Compute pressure from D3Q19 distribution
* @param Dist - D3Q19 distributions
* \brief compute pressure from D3Q19 distribution
* @param dist - D3Q19 distributions
* @param press - memory buffer to store the pressure field that is computed
* @param Np - size of local sub-domain (derived from Domain structure)
*/
@ -250,12 +250,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);
@ -269,10 +269,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){

15
docs/README.md
View File

@ -4,12 +4,15 @@ Dependencies for LBPM documentation
```python
pip install Sphinx
```
# foamatting requires sphinx read-the-docs-theme
# formatting requires sphinx read-the-docs-theme
```python
pip install sphinx-rtd-theme
```
# support for doxygen requires breathe
pip install breathe
# equation rendering requires latex and dvipng command
```
sudo apt-get install dvipng
@ -17,11 +20,13 @@ sudo apt-get install texlive texstudio
sudo apt-get install texlive-latex-recommended texlive-pictures texlive-latex-extra
```
# To build the docs
```
Step 1) install dependencies listed above
Step 2) type 'make html' from the docs/ directory
Step 3) point your browser at ~/local/doc/build/html/index.html
```
Step 2) build LBPM with doxygen support enabled in CMake
Step 3) From the build directory run the command 'make doc' to build html and xml from doxygen
Step 4) modify source/conf.py so that the breathe directory specifies the full path to the xml
documentation built in Step 3
Step 5) type 'make html' from the docs/ directory (the same directory that contains this README file)
Step 6) point your browser at ~/local/doc/build/html/index.html
#

13
docs/source/conf.py
View File

@ -31,7 +31,8 @@ release = '1.0'
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones.
extensions = [
'sphinx.ext.imgmath'
'sphinx.ext.imgmath',
'breathe'
]
# Add any paths that contain templates here, relative to this directory.
@ -55,7 +56,6 @@ html_theme = 'alabaster'
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
## Read the docs style:
if os.environ.get('READTHEDOCS') != 'True':
try:
@ -68,3 +68,12 @@ if os.environ.get('READTHEDOCS') != 'True':
#def setup(app):
# app.add_stylesheet("fix_rtd.css")
# -- Breathe configuration -------------------------------------------------
breathe_projects = {
"LBPM Doxygen": "/home/mcclurej/local/dev/LBPM/doc/xml/"
}
breathe_default_project = "LBPM Doxygen"
breathe_default_members = ('members', 'undoc-members')

9
docs/source/developerGuide/doxygen/Analysis.rst Normal file
View File

@ -0,0 +1,9 @@
###############################################################################
Analysis methods
###############################################################################
.. toctree::
:glob:
:maxdepth: 2
analysis/*

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

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

10
docs/source/developerGuide/doxygen/Models.rst Normal file
View File

@ -0,0 +1,10 @@
###############################################################################
Models
###############################################################################
.. toctree::
:glob:
:maxdepth: 2
models/*

5
docs/source/developerGuide/doxygen/ScaLBL.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
Scalable Lattice Boltzmann Library (ScalBL)
============================================
.. doxygenfile:: ScaLBL.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/analysis/DCEL.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
DCEL
============================================
.. doxygenfile:: dcel.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/analysis/FlowAdaptor.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
FlowAdaptor
============================================
.. doxygenfile:: FlowAdaptor.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/analysis/Minkowski.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
Minkowski
============================================
.. doxygenfile:: Minkowski.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/analysis/SubPhase.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
SubPhase analysis
============================================
.. doxygenfile:: SubPhase.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/models/ColorModel.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
Color Model
============================================
.. doxygenfile:: ColorModel.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/models/FreeEnergyModel.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
Free Energy Model
============================================
.. doxygenfile:: FreeLeeModel.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/models/GreyscaleModel.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
Greyscale Model
============================================
.. doxygenfile:: GreyscaleModel.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/models/MRTModel.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
Multi-relaxation time (MRT) Model
============================================
.. doxygenfile:: MRTModel.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/models/NernstPlanckModel.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
Ion Model
============================================
.. doxygenfile:: IonModel.h
:project: LBPM Doxygen

5
docs/source/developerGuide/doxygen/models/PoissonModel.rst Normal file
View File

@ -0,0 +1,5 @@
============================================
Poisson Model
============================================
.. doxygenfile:: PoissonSolver.h
:project: LBPM Doxygen

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

@ -15,4 +15,10 @@ into the framework.
testingModels/*
doxygen/ScaLBL.rst
doxygen/Domain.rst
doxygen/Models.rst
doxygen/Analysis.rst

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:`\epsilon_m`, 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/*

2
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`

6
doxygen/Doxyfile.in
View File

@ -25,13 +25,13 @@ DOXYFILE_ENCODING = UTF-8
# The PROJECT_NAME tag is a single word (or a sequence of words surrounded
# by quotes) that should identify the project.
PROJECT_NAME = "@PROJECT_NAME@"
PROJECT_NAME = "LBPM Doxygen"
# The PROJECT_NUMBER tag can be used to enter a project or revision number.
# This could be handy for archiving the generated documentation or
# if some version control system is used.
PROJECT_NUMBER = "@PROJECT_VERSION@"
PROJECT_NUMBER = "9.12.2021"
# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewer
@ -985,7 +985,7 @@ MAN_LINKS = NO
# generate an XML file that captures the structure of
# the code including all documentation.
GENERATE_XML = NO
GENERATE_XML = YES
# The XML_OUTPUT tag is used to specify where the XML pages will be put.
# If a relative path is entered the value of OUTPUT_DIRECTORY will be

3
doxygen/DoxygenMainpage.h
View File

@ -3,8 +3,9 @@
* C/C++ routines
*
* - \ref ScaLBL.h "Scalable Lattice Boltzmann Library (ScaLBL)"
* - \ref Domain.h "Domain structure"
* - \ref models "Lattice Boltzmann models"
* - \ref ColorModel "Color model"
* - \ref ScaLBL_ColorModel "Color model"
* - \ref analysis "Analysis routines"
* - \ref FlowAdaptor "FlowAdaptor"
* - \ref DCEL "Doubly connected edge list"

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

@ -127,7 +127,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;
//}
//----------------------------------------------------------------------------------------------------------//

70
models/GreyscaleColorModel.h
View File

@ -15,19 +15,69 @@ 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();
bool Restart,pBC;
@ -72,7 +122,7 @@ public:
double *GreySw;
double *GreyKn;
double *GreyKw;
//double *ColorGrad;
double *MobilityRatio;
double *Velocity;
double *Pressure;
double *Porosity_dvc;
@ -91,14 +141,24 @@ 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

@ -632,29 +632,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++){
@ -781,51 +758,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.WriteVis();
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_greyscaleColor_simulator",1);