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

merge master into FOM_dev

Browse Source
This commit is contained in:
Rex Zhe Li 2021-06-27 19:53:52 -04:00
commit de909820f8
58 changed files with 2650 additions and 791 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists.txt
View File

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

49
analysis/ElectroChemistry.cpp
View File

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

3
analysis/FreeEnergy.cpp
View File

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

196
analysis/SubPhase.cpp
View File

@ -23,6 +23,7 @@ SubPhase::SubPhase(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);
Dissipation.resize(Nx,Ny,Nz); Dissipation.fill(0);
SDs.resize(Nx,Ny,Nz); SDs.fill(0);
//.........................................
@ -42,11 +43,12 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");
fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn wet ");
fprintf(SUBPHASE,"pwc pwd pnc pnd "); // pressures
fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni "); // mass
fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x "); // momentum
fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");
fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");
fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni Msw Msn "); // mass
fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x Psw_x Psn_x "); // momentum
fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y Psw_y Psn_y ");
fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z Psw_z Psn_z ");
fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy
fprintf(SUBPHASE,"Dwc Dwd Dnc Dnd "); // viscous dissipation
fprintf(SUBPHASE,"Vwc Awc Hwc Xwc "); // wc region
fprintf(SUBPHASE,"Vwd Awd Hwd Xwd Nwd "); // wd region
fprintf(SUBPHASE,"Vnc Anc Hnc Xnc "); // nc region
@ -67,11 +69,12 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");
fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn wet ");
fprintf(SUBPHASE,"pwc pwd pnc pnd "); // pressures
fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni "); // mass
fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x "); // momentum
fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");
fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");
fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni Msw Msn "); // mass
fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x Psw_x Psn_x "); // momentum
fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y Psw_y Psn_y ");
fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z Psw_z Psn_z ");
fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy
fprintf(SUBPHASE,"Dwc Dwd Dnc Dnd "); // viscous dissipation
fprintf(SUBPHASE,"Vwc Awc Hwc Xwc "); // wc region
fprintf(SUBPHASE,"Vwd Awd Hwd Xwd Nwd "); // wd region
fprintf(SUBPHASE,"Vnc Anc Hnc Xnc "); // nc region
@ -93,7 +96,7 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
{
// If timelog is empty, write a short header to list the averages
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"sw krw krn vw vn force pw pn wet\n");
fprintf(TIMELOG,"sw krw krn krwf krnf vw vn force pw pn wet\n");
}
}
}
@ -111,11 +114,12 @@ void SubPhase::Write(int timestep)
if (Dm->rank()==0){
fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn,total_wetting_interaction_global);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.p, gwd.p, gnc.p, gnd.p);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn, gifs.Mw, gifs.Mn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx, gifs.Pwx, gifs.Pnx);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny, gifs.Pwy, gifs.Pny);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz, gifs.Pwz, gifs.Pnz);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.visc, gwd.visc, gnc.visc, gnd.visc);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.V, gwc.A, gwc.H, gwc.X);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gnc.V, gnc.A, gnc.H, gnc.X);
@ -127,11 +131,12 @@ void SubPhase::Write(int timestep)
else{
fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn,total_wetting_interaction);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.p, wd.p, nc.p, nd.p);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn, ifs.Mw, ifs.Mn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx, ifs.Pwx, ifs.Pnx);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny, ifs.Pwy, ifs.Pny);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz, ifs.Pwz, ifs.Pnz);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.visc, wd.visc, nc.visc, nd.visc);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.V, wc.A, wc.H, wc.X);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",nc.V, nc.A, nc.H, nc.X);
@ -167,7 +172,7 @@ void SubPhase::Basic(){
if (Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin += Dm->inlet_layers_z;
if (Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax -= Dm->outlet_layers_z;
*/
nb.reset(); wb.reset();
nb.reset(); wb.reset(); iwn.reset();
double count_w = 0.0;
double count_n = 0.0;
@ -245,6 +250,27 @@ void SubPhase::Basic(){
wb.p += Pressure(n);
count_w += 1.0;
}
/* compute the film contribution */
else if (SDs(i,j,k) < 2.0){
if ( phi > 0.0 ){
nA = 1.0;
iwn.V += 1.0;
iwn.Mn += nA*rho_n;
// velocity
iwn.Pnx += rho_n*nA*Vel_x(n);
iwn.Pny += rho_n*nA*Vel_y(n);
iwn.Pnz += rho_n*nA*Vel_z(n);
}
else{
nB = 1.0;
iwn.Mw += nB*rho_w;
iwn.V += 1.0;
iwn.Pwx += rho_w*nB*Vel_x(n);
iwn.Pwy += rho_w*nB*Vel_y(n);
iwn.Pwz += rho_w*nB*Vel_z(n);
}
}
}
}
}
@ -267,12 +293,6 @@ void SubPhase::Basic(){
//printf("wetting interaction = %f, count = %f\n",total_wetting_interaction,count_wetting_interaction);
total_wetting_interaction_global=Dm->Comm.sumReduce( total_wetting_interaction);
count_wetting_interaction_global=Dm->Comm.sumReduce( count_wetting_interaction);
/* normalize wetting interactions <-- Don't do this if normalizing laplacian (use solid surface area)
if (count_wetting_interaction > 0.0)
total_wetting_interaction /= count_wetting_interaction;
if (count_wetting_interaction_global > 0.0)
total_wetting_interaction_global /= count_wetting_interaction_global;
*/
gwb.V=Dm->Comm.sumReduce( wb.V);
gnb.V=Dm->Comm.sumReduce( nb.V);
@ -285,6 +305,16 @@ void SubPhase::Basic(){
gnb.Py=Dm->Comm.sumReduce( nb.Py);
gnb.Pz=Dm->Comm.sumReduce( nb.Pz);
giwn.Mw=Dm->Comm.sumReduce( iwn.Mw);
giwn.Pwx=Dm->Comm.sumReduce( iwn.Pwx);
giwn.Pwy=Dm->Comm.sumReduce( iwn.Pwy);
giwn.Pwz=Dm->Comm.sumReduce( iwn.Pwz);
giwn.Mn=Dm->Comm.sumReduce( iwn.Mn);
giwn.Pnx=Dm->Comm.sumReduce( iwn.Pnx);
giwn.Pny=Dm->Comm.sumReduce( iwn.Pny);
giwn.Pnz=Dm->Comm.sumReduce( iwn.Pnz);
count_w=Dm->Comm.sumReduce( count_w);
count_n=Dm->Comm.sumReduce( count_n);
if (count_w > 0.0)
@ -310,20 +340,21 @@ void SubPhase::Basic(){
if (gnb.Pz != gnb.Pz) err=true;
if (Dm->rank() == 0){
/* align flow direction based on total mass flux */
double dir_x = gwb.Px + gnb.Px;
double dir_y = gwb.Py + gnb.Py;
double dir_z = gwb.Pz + gnb.Pz;
double flow_magnitude = dir_x*dir_x + dir_y*dir_y + dir_z*dir_z;
double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
double dir_x = 0.0;
double dir_y = 0.0;
double dir_z = 0.0;
if (force_mag > 0.0){
dir_x = Fx/force_mag;
dir_y = Fy/force_mag;
dir_z = Fz/force_mag;
}
else {
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
dir_z = 1.0;
dir_x /= flow_magnitude;
dir_y /= flow_magnitude;
dir_z /= flow_magnitude;
}
if (Dm->BoundaryCondition == 1 || Dm->BoundaryCondition == 2 || Dm->BoundaryCondition == 3 || Dm->BoundaryCondition == 4 ){
// compute the pressure drop
@ -331,7 +362,7 @@ void SubPhase::Basic(){
double length = ((Nz-2)*Dm->nprocz());
force_mag -= pressure_drop/length;
}
if (force_mag == 0.0){
if (force_mag == 0.0 && flow_magnitude == 0.0){
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
@ -341,14 +372,21 @@ void SubPhase::Basic(){
double saturation=gwb.V/(gwb.V + gnb.V);
double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume;
double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume;
/* contribution from water films */
double water_film_flow_rate=gwb.V*(giwn.Pwx*dir_x + giwn.Pwy*dir_y + giwn.Pwz*dir_z)/gwb.M / Dm->Volume;
double not_water_film_flow_rate=gnb.V*(giwn.Pnx*dir_x + giwn.Pny*dir_y + giwn.Pnz*dir_z)/gnb.M / Dm->Volume;
//double total_flow_rate = water_flow_rate + not_water_flow_rate;
//double fractional_flow = water_flow_rate / total_flow_rate;
double h = Dm->voxel_length;
double krn = h*h*nu_n*not_water_flow_rate / force_mag ;
double krw = h*h*nu_w*water_flow_rate / force_mag;
/* not counting films */
double krnf = krn - h*h*nu_n*not_water_film_flow_rate / force_mag ;
double krwf = krw - h*h*nu_w*water_film_flow_rate / force_mag;
//printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);
fprintf(TIMELOG,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",saturation,krw,krn,h*water_flow_rate,h*not_water_flow_rate, force_mag, gwb.p, gnb.p, total_wetting_interaction_global);
fprintf(TIMELOG,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",saturation,krw,krn,krwf,krnf,h*water_flow_rate,h*not_water_flow_rate, force_mag, gwb.p, gnb.p, total_wetting_interaction_global);
fflush(TIMELOG);
}
if (err==true){
@ -364,6 +402,7 @@ inline void InterfaceTransportMeasures( double beta, double rA, double rB, doubl
double A1,A2,A3,A4,A5,A6;
double B1,B2,B3,B4,B5,B6;
double nAB,delta;
double phi = (nA-nB)/(nA+nB);
// Instantiate mass transport distributions
// Stationary value - distribution 0
nAB = 1.0/(nA+nB);
@ -402,7 +441,7 @@ inline void InterfaceTransportMeasures( double beta, double rA, double rB, doubl
double uwx = (B1-B2);
double uwy = (B3-B4);
double uwz = (B5-B6);
/*
I.Mn += rA*nA;
I.Mw += rB*nB;
I.Pnx += rA*nA*unx;
@ -413,6 +452,21 @@ inline void InterfaceTransportMeasures( double beta, double rA, double rB, doubl
I.Pwz += rB*nB*uwz;
I.Kn += rA*nA*(unx*unx + uny*uny + unz*unz);
I.Kw += rB*nB*(uwx*uwx + uwy*uwy + uwz*uwz);
*/
if (phi > 0.0){
I.Mn += rA;
I.Pnx += rA*ux;
I.Pny += rA*uy;
I.Pnz += rA*uz;
}
else {
I.Mw += rB;
I.Pwx += rB*ux;
I.Pwy += rB*uy;
I.Pwz += rB*uz;
}
I.Kn += rA*nA*(unx*unx + uny*uny + unz*unz);
I.Kw += rB*nB*(uwx*uwx + uwy*uwy + uwz*uwz);
}
@ -432,7 +486,7 @@ void SubPhase::Full(){
if (Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin += Dm->inlet_layers_z;
if (Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax -= Dm->outlet_layers_z;
*/
nd.reset(); nc.reset(); wd.reset(); wc.reset(); iwn.reset(); iwnc.reset();
nd.reset(); nc.reset(); wd.reset(); wc.reset(); iwn.reset(); iwnc.reset(); ifs.reset();
Dm->CommunicateMeshHalo(Phi);
for (int k=1; k<Nz-1; k++){
@ -447,6 +501,33 @@ void SubPhase::Full(){
}
}
Dm->CommunicateMeshHalo(DelPhi);
Dm->CommunicateMeshHalo(Vel_x);
Dm->CommunicateMeshHalo(Vel_y);
Dm->CommunicateMeshHalo(Vel_z);
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
// Compute velocity gradients using finite differences
double phi = Phi(i,j,k);
double nu = nu_n + 0.5*(1.0-phi)*(nu_w-nu_n);
double rho = rho_n + 0.5*(1.0-phi)*(rho_w-rho_n);
double ux = 0.5*(Vel_x(i+1,j,k) - Vel_x(i-1,j,k));
double uy = 0.5*(Vel_x(i,j+1,k) - Vel_x(i,j-1,k));
double uz = 0.5*(Vel_x(i,j,k+1) - Vel_x(i,j,k-1));
double vx = 0.5*(Vel_y(i+1,j,k) - Vel_y(i-1,j,k));
double vy = 0.5*(Vel_y(i,j+1,k) - Vel_y(i,j-1,k));
double vz = 0.5*(Vel_y(i,j,k+1) - Vel_y(i,j,k-1));
double wx = 0.5*(Vel_z(i+1,j,k) - Vel_z(i-1,j,k));
double wy = 0.5*(Vel_z(i,j+1,k) - Vel_z(i,j-1,k));
double wz = 0.5*(Vel_z(i,j,k+1) - Vel_z(i,j,k-1));
if (SDs(i,j,k) > 2.0){
Dissipation(i,j,k) = 2*rho*nu*( ux*ux + vy*vy + wz*wz + 0.5*(vx + uy)*(vx + uy)+ 0.5*(vz + wy)*(vz + wy)+ 0.5*(uz + wx)*(uz + wx));
}
}
}
}
/* Set up geometric analysis of each region */
@ -605,13 +686,33 @@ void SubPhase::Full(){
double ux = Vel_x(n);
double uy = Vel_y(n);
double uz = Vel_z(n);
double visc = Dissipation(n);
if (DelPhi(n) > 1e-3){
// interface region
if (DelPhi(n) > 1e-3 ){
// get the normal vector
double nx = 0.5*(Phi(i+1,j,k)-Phi(i-1,j,k));
double ny = 0.5*(Phi(i,j+1,k)-Phi(i,j-1,k));
double nz = 0.5*(Phi(i,j,k+1)-Phi(i,j,k-1));
InterfaceTransportMeasures( beta, rho_w, rho_n, nA, nB, nx, ny, nz, ux, uy, uz, iwn);
if (SDs(n) > 2.5){
// not a film region
InterfaceTransportMeasures( beta, rho_w, rho_n, nA, nB, nx, ny, nz, ux, uy, uz, iwn);
}
else{
// films that are close to the wetting fluid
if ( morph_w->distance(i,j,k) < 2.5 && phi > 0.0){
ifs.Mw += rho_w;
ifs.Pwx += rho_w*ux;
ifs.Pwy += rho_w*uy;
ifs.Pwz += rho_w*uz;
}
// films that are close to the NWP
if ( morph_n->distance(i,j,k) < 2.5 && phi < 0.0){
ifs.Mn += rho_n;
ifs.Pnx += rho_n*ux;
ifs.Pny += rho_n*uy;
ifs.Pnz += rho_n*uz;
}
}
}
else if ( phi > 0.0){
if (morph_n->label(i,j,k) > 0 ){
@ -642,6 +743,7 @@ void SubPhase::Full(){
nd.Py += nA*rho_n*uy;
nd.Pz += nA*rho_n*uz;
nd.K += nA*rho_n*(ux*ux + uy*uy + uz*uz);
nd.visc += visc;
}
else{
nA = 1.0;
@ -650,6 +752,7 @@ void SubPhase::Full(){
nc.Py += nA*rho_n*uy;
nc.Pz += nA*rho_n*uz;
nc.K += nA*rho_n*(ux*ux + uy*uy + uz*uz);
nc.visc += visc;
}
}
else{
@ -661,6 +764,7 @@ void SubPhase::Full(){
wd.Py += nB*rho_w*uy;
wd.Pz += nB*rho_w*uz;
wd.K += nB*rho_w*(ux*ux + uy*uy + uz*uz);
wd.visc += visc;
}
else{
nB = 1.0;
@ -669,6 +773,7 @@ void SubPhase::Full(){
wc.Py += nB*rho_w*uy;
wc.Pz += nB*rho_w*uz;
wc.K += nB*rho_w*(ux*ux + uy*uy + uz*uz);
wc.visc += visc;
}
}
}
@ -681,25 +786,29 @@ void SubPhase::Full(){
gnd.Py=Dm->Comm.sumReduce( nd.Py);
gnd.Pz=Dm->Comm.sumReduce( nd.Pz);
gnd.K=Dm->Comm.sumReduce( nd.K);
gnd.visc=Dm->Comm.sumReduce( nd.visc);
gwd.M=Dm->Comm.sumReduce( wd.M);
gwd.Px=Dm->Comm.sumReduce( wd.Px);
gwd.Py=Dm->Comm.sumReduce( wd.Py);
gwd.Pz=Dm->Comm.sumReduce( wd.Pz);
gwd.K=Dm->Comm.sumReduce( wd.K);
gwd.visc=Dm->Comm.sumReduce( wd.visc);
gnc.M=Dm->Comm.sumReduce( nc.M);
gnc.Px=Dm->Comm.sumReduce( nc.Px);
gnc.Py=Dm->Comm.sumReduce( nc.Py);
gnc.Pz=Dm->Comm.sumReduce( nc.Pz);
gnc.K=Dm->Comm.sumReduce( nc.K);
gnc.visc=Dm->Comm.sumReduce( nc.visc);
gwc.M=Dm->Comm.sumReduce( wc.M);
gwc.Px=Dm->Comm.sumReduce( wc.Px);
gwc.Py=Dm->Comm.sumReduce( wc.Py);
gwc.Pz=Dm->Comm.sumReduce( wc.Pz);
gwc.K=Dm->Comm.sumReduce( wc.K);
gwc.visc=Dm->Comm.sumReduce( wc.visc);
giwn.Mn=Dm->Comm.sumReduce( iwn.Mn);
giwn.Pnx=Dm->Comm.sumReduce( iwn.Pnx);
giwn.Pny=Dm->Comm.sumReduce( iwn.Pny);
@ -711,6 +820,15 @@ void SubPhase::Full(){
giwn.Pwz=Dm->Comm.sumReduce( iwn.Pwz);
giwn.Kw=Dm->Comm.sumReduce( iwn.Kw);
gifs.Mn= Dm->Comm.sumReduce( ifs.Mn);
gifs.Pnx=Dm->Comm.sumReduce( ifs.Pnx);
gifs.Pny=Dm->Comm.sumReduce( ifs.Pny);
gifs.Pnz=Dm->Comm.sumReduce( ifs.Pnz);
gifs.Mw= Dm->Comm.sumReduce( ifs.Mw);
gifs.Pwx=Dm->Comm.sumReduce( ifs.Pwx);
gifs.Pwy=Dm->Comm.sumReduce( ifs.Pwy);
gifs.Pwz=Dm->Comm.sumReduce( ifs.Pwz);
// pressure averaging
gnc.p=Dm->Comm.sumReduce( nc.p);
gnd.p=Dm->Comm.sumReduce( nd.p);

6
analysis/SubPhase.h
View File

@ -22,10 +22,11 @@ class phase{
public:
int Nc;
double p;
double M,Px,Py,Pz,K;
double M,Px,Py,Pz,K,visc;
double V,A,H,X;
void reset(){
p=M=Px=Py=Pz=K=0.0;
visc=0.0;
V=A=H=X=0.0;
Nc=1;
}
@ -70,10 +71,12 @@ public:
// local entities
phase wc,wd,wb,nc,nd,nb,solid;
interface iwn,iwnc;
interface ifs;
// global entities
phase gwc,gwd,gwb,gnc,gnd,gnb,gsolid;
interface giwn,giwnc;
interface gifs;
/* fluid-solid wetting interaction */
double total_wetting_interaction, count_wetting_interaction;
double total_wetting_interaction_global, count_wetting_interaction_global;
@ -92,6 +95,7 @@ public:
DoubleArray Vel_x; // velocity field
DoubleArray Vel_y;
DoubleArray Vel_z;
DoubleArray Dissipation;
DoubleArray SDs;
std::shared_ptr<Minkowski> morph_w;

20
analysis/morphology.cpp
View File

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

24
analysis/runAnalysis.cpp
View File

@ -296,15 +296,21 @@ public:
fillData.copy( Averages.Vel_z, VelzData );
}
if ( vis_db->getWithDefault<bool>( "save_dissipation", false ) ) {
ASSERT( visData[0].vars[5]->name == "ViscousDissipation" );
Array<double> &ViscousDissipation = visData[0].vars[5]->data;
fillData.copy( Averages.Dissipation, ViscousDissipation );
}
if ( vis_db->getWithDefault<bool>( "save_distance", false ) ) {
ASSERT( visData[0].vars[5]->name == "SignDist" );
Array<double> &SignData = visData[0].vars[5]->data;
ASSERT( visData[0].vars[6]->name == "SignDist" );
Array<double> &SignData = visData[0].vars[6]->data;
fillData.copy( Averages.SDs, SignData );
}
if ( vis_db->getWithDefault<bool>( "save_connected_components", false ) ) {
ASSERT( visData[0].vars[6]->name == "BlobID" );
Array<double> &BlobData = visData[0].vars[6]->data;
ASSERT( visData[0].vars[7]->name == "BlobID" );
Array<double> &BlobData = visData[0].vars[7]->data;
fillData.copy( Averages.morph_n->label, BlobData );
}
@ -652,6 +658,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
auto ViscousDissipationVar = std::make_shared<IO::Variable>();
auto SignDistVar = std::make_shared<IO::Variable>();
auto BlobIDVar = std::make_shared<IO::Variable>();
@ -688,6 +695,14 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
VzVar->data.resize( d_n[0], d_n[1], d_n[2] );
d_meshData[0].vars.push_back( VzVar );
}
if ( vis_db->getWithDefault<bool>( "save_dissipation", false ) ) {
ViscousDissipationVar->name = "ViscousDissipation";
ViscousDissipationVar->type = IO::VariableType::VolumeVariable;
ViscousDissipationVar->dim = 1;
ViscousDissipationVar->data.resize( d_n[0], d_n[1], d_n[2] );
d_meshData[0].vars.push_back( ViscousDissipationVar );
}
if ( vis_db->getWithDefault<bool>( "save_distance", false ) ) {
SignDistVar->name = "SignDist";
@ -729,7 +744,6 @@ runAnalysis::runAnalysis( ScaLBL_ColorModel &ColorModel)
d_comm = ColorModel.Dm->Comm.dup();
d_Np = ColorModel.Np;
bool Regular = false;
auto input_db = ColorModel.db;
auto db = input_db->getDatabase( "Analysis" );

2
common/Domain.cpp
View File

@ -383,7 +383,7 @@ void Domain::Decomp( const std::string& Filename )
for (int i = 0; i<global_Nx; i++){
n = k*global_Nx*global_Ny+j*global_Nx+i;
//char locval = loc_id[n];
char locval = SegData[n];
signed char locval = SegData[n];
for (size_t idx=0; idx<ReadValues.size(); idx++){
signed char oldvalue=ReadValues[idx];
signed char newvalue=WriteValues[idx];

135
common/ScaLBL.cpp
View File

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

4
common/ScaLBL.h
View File

@ -88,12 +88,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 *Poros,double *Perm,double *Vel, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, 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 *Poros,double *Perm,double *Vel,double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);

5
cpu/BGK.cpp
View File

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

20
cpu/Color.cpp
View File

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

12
cpu/Greyscale.cpp
View File

@ -2,7 +2,6 @@
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity, double *Pressure){
int n;
// conserved momemnts
double rho,vx,vy,vz,v_mag;
double ux,uy,uz,u_mag;
@ -247,7 +246,6 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finis
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity,double *Pressure){
int n;
// conserved momemnts
double rho,vx,vy,vz,v_mag;
double ux,uy,uz,u_mag;
@ -263,7 +261,6 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, in
double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
int nread;
for (int n=start; n<finish; n++){
// q=0
@ -553,7 +550,6 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, in
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity, double Den,double *Pressure){
int n;
double vx,vy,vz,v_mag;
double ux,uy,uz,u_mag;
double pressure;//defined for this incompressible model
@ -1042,7 +1038,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,
double *Poros,double *Perm, double *Velocity, double Den,double *Pressure){
int n, nread;
int nread;
double vx,vy,vz,v_mag;
double ux,uy,uz,u_mag;
double pressure;//defined for this incompressible model
@ -1197,6 +1193,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
// q=9
nread = neighborList[n+8*Np];
fq = dist[nread];
pressure += fq;
m1 += 8.0*fq;
@ -1568,7 +1565,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,double *Poros,double *Perm, double *Velocity,double rho0,double *Pressure){
int n, nread;
int nread;
int nr1,nr2,nr3,nr4,nr5,nr6;
int nr7,nr8,nr9,nr10;
int nr11,nr12,nr13,nr14;
@ -1705,6 +1702,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
//nread = neighborList[n+6*Np];
//fq = dist[nread];
nr7 = neighborList[n+6*Np];
fq = dist[nr7];
rho += fq;
m1 += 8.0*fq;
@ -1954,6 +1952,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
if (porosity==1.0){
Fx=rho0*Gx;
Fy=rho0*Gy;
Fz=rho0*Gz;
}
@ -2120,7 +2119,6 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist
extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_MRT(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,double *Poros,double *Perm, double *Velocity,double rho0,double *Pressure){
int n;
double vx,vy,vz,v_mag;
double ux,uy,uz,u_mag;
double pressure;//defined for this incompressible model

263
cpu/GreyscaleColor.cpp
View File

@ -1340,10 +1340,10 @@ 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 *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Velocity,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, double W, int strideY, int strideZ, int start, int finish, int Np){
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 *Poros,double *Perm, double *Velocity, 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){
int n,nn,ijk,nread;
int nr1,nr2,nr3,nr4,nr5,nr6;
@ -1370,12 +1370,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//double c0, c1; //Guo's model parameters
double tau_eff;
double mu_eff;//kinematic viscosity
double nx_gs,ny_gs,nz_gs;//grey-solid color gradient
double nx_phase,ny_phase,nz_phase,C_phase;
double Fx,Fy,Fz;
double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18;
double gp3,gp5,gp7;
double Fcpx,Fcpy,Fcpz;//capillary penalty force
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
const double mrt_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802;
@ -1394,11 +1392,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
// read the component number densities
nA = Den[n];
nB = Den[Np + n];
porosity = Poros[n];
perm = Perm[n];
nx_gs = GreySolidGrad[n+0*Np];
ny_gs = GreySolidGrad[n+1*Np];
nz_gs = GreySolidGrad[n+2*Np];
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
@ -1420,100 +1419,63 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
gp1 = Psi[nn];
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
gp2 = Psi[nn];
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Phi[nn]; // get neighbor for phi - 3
gp3 = Psi[nn];
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Phi[nn]; // get neighbor for phi - 4
gp4 = Psi[nn];
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Phi[nn]; // get neighbor for phi - 5
gp5 = Psi[nn];
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Phi[nn]; // get neighbor for phi - 6
gp6 = Psi[nn];
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Phi[nn]; // get neighbor for phi - 7
gp7 = Psi[nn];
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Phi[nn]; // get neighbor for phi - 8
gp8 = Psi[nn];
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Phi[nn]; // get neighbor for phi - 9
gp9 = Psi[nn];
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Phi[nn]; // get neighbor for phi - 10
gp10 = Psi[nn];
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Phi[nn]; // get neighbor for phi - 11
gp11 = Psi[nn];
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Phi[nn]; // get neighbor for phi - 12
gp12 = Psi[nn];
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Phi[nn]; // get neighbor for phi - 13
gp13 = Psi[nn];
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Phi[nn]; // get neighbor for phi - 14
gp14 = Psi[nn];
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Phi[nn]; // get neighbor for phi - 15
gp15 = Psi[nn];
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Phi[nn]; // get neighbor for phi - 16
gp16 = Psi[nn];
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Phi[nn]; // get neighbor for phi - 17
gp17 = Psi[nn];
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Phi[nn]; // get neighbor for phi - 18
gp18 = Psi[nn];
//............Compute the Color Gradient...................................
nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase);
nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//correct the normal color gradient by considering the effect of grey solid
nx = nx_phase + (1.0-porosity)*nx_gs;
ny = ny_phase + (1.0-porosity)*ny_gs;
nz = nz_phase + (1.0-porosity)*nz_gs;
if (C_phase==0.0){
nx = nx_phase;
ny = ny_phase;
nz = nz_phase;
}
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//............Compute the Greyscale Potential Gradient.....................
//............Compute the Greyscale Potential Gradient.....................
// Fcpx = 0.0;
// Fcpy = 0.0;
// Fcpz = 0.0;
@ -1534,14 +1496,24 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
// //ny = Fcpy/Fcp_mag;
// //nz = Fcpz/Fcp_mag;
// }
Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
Fcpx = nx;
Fcpy = ny;
Fcpz = nz;
double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
if (Fcp_mag==0.0) Fcpx=Fcpy=Fcpz=0.0;
//NOTE for open node (porosity=1.0),Fcp=0.0
Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm);
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=0
fq = dist[n];
rho = fq;
@ -1893,6 +1865,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//........................................................................
//..............carry out relaxation process..............................
//..........Toelke, Fruediger et. al. 2006................................
//---------------- NO higher-order force -------------------------------//
if (C == 0.0) nx = ny = nz = 0.0;
m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2);
@ -1917,6 +1890,43 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//----------------With higher-order force ------------------------------//
//if (C == 0.0) nx = ny = nz = 0.0;
//m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1)
// + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
//m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2)
// + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity;
//jx = jx + Fx;
//m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
//jy = jy + Fy;
//m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
//jz = jz + Fz;
//m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
//m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9)
// + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
////m10 = m10 + rlx_setA*( - m10);
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10)
// + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity;
//m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11)
// + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
////m12 = m12 + rlx_setA*( - m12);
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12)
// + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity;
//m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13);
// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
//m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14);
// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
//m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15);
// + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity;
//m16 = m16 + rlx_setB*( - m16);
//m17 = m17 + rlx_setB*( - m17);
//m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//.................inverse transformation......................................................
// q=0
@ -2049,6 +2059,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@ -2068,6 +2082,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
// Cq = {0,1,0}
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@ -2088,6 +2106,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
// Cq = {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@ -2109,9 +2131,9 @@ 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 *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Velocity,double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, 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, double W, int strideY, int strideZ, int start, int finish, int Np){
double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
int ijk,nn,n;
double fq;
@ -2127,18 +2149,16 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
double a1,b1,a2,b2,nAB,delta;
double C,nx,ny,nz; //color gradient magnitude and direction
double phi,tau,rho0,rlx_setA,rlx_setB;
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
//double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002)
double porosity;
double perm;//voxel permeability
//double c0, c1; //Guo's model parameters
double tau_eff;
double mu_eff;//kinematic viscosity
double nx_gs,ny_gs,nz_gs;//grey-solid color gradient
double nx_phase,ny_phase,nz_phase,C_phase;
double Fx,Fy,Fz;
double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18;
double gp3,gp5,gp7;
double Fcpx,Fcpy,Fcpz;//capillary penalty force
const double mrt_V1=0.05263157894736842;
@ -2155,15 +2175,15 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
const double mrt_V12=0.04166666666666666;
for (n=start; n<finish; n++){
// read the component number densities
nA = Den[n];
nB = Den[Np + n];
porosity = Poros[n];
perm = Perm[n];
nx_gs = GreySolidGrad[n+0*Np];
ny_gs = GreySolidGrad[n+1*Np];
nz_gs = GreySolidGrad[n+2*Np];
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
@ -2185,100 +2205,63 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
gp1 = Psi[nn];
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
gp2 = Psi[nn];
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Phi[nn]; // get neighbor for phi - 3
gp3 = Psi[nn];
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Phi[nn]; // get neighbor for phi - 4
gp4 = Psi[nn];
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Phi[nn]; // get neighbor for phi - 5
gp5 = Psi[nn];
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Phi[nn]; // get neighbor for phi - 6
gp6 = Psi[nn];
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Phi[nn]; // get neighbor for phi - 7
gp7 = Psi[nn];
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Phi[nn]; // get neighbor for phi - 8
gp8 = Psi[nn];
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Phi[nn]; // get neighbor for phi - 9
gp9 = Psi[nn];
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Phi[nn]; // get neighbor for phi - 10
gp10 = Psi[nn];
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Phi[nn]; // get neighbor for phi - 11
gp11 = Psi[nn];
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Phi[nn]; // get neighbor for phi - 12
gp12 = Psi[nn];
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Phi[nn]; // get neighbor for phi - 13
gp13 = Psi[nn];
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Phi[nn]; // get neighbor for phi - 14
gp14 = Psi[nn];
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Phi[nn]; // get neighbor for phi - 15
gp15 = Psi[nn];
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Phi[nn]; // get neighbor for phi - 16
gp16 = Psi[nn];
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Phi[nn]; // get neighbor for phi - 17
gp17 = Psi[nn];
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Phi[nn]; // get neighbor for phi - 18
gp18 = Psi[nn];
//............Compute the Color Gradient...................................
nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase);
nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//correct the normal color gradient by considering the effect of grey solid
nx = nx_phase + (1.0-porosity)*nx_gs;
ny = ny_phase + (1.0-porosity)*ny_gs;
nz = nz_phase + (1.0-porosity)*nz_gs;
if (C_phase==0.0){
nx = nx_phase;
ny = ny_phase;
nz = nz_phase;
}
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
//............Compute the Greyscale Potential Gradient.....................
//............Compute the Greyscale Potential Gradient.....................
// Fcpx = 0.0;
// Fcpy = 0.0;
// Fcpz = 0.0;
@ -2299,14 +2282,23 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
// ny = Fcpy/Fcp_mag;
// nz = Fcpz/Fcp_mag;
// }
Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
Fcpx = nx;
Fcpy = ny;
Fcpz = nz;
double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
if (Fcp_mag==0.0) Fcpx=Fcpy=Fcpz=0.0;
//NOTE for open node (porosity=1.0),Fcp=0.0
Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm);
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=0
fq = dist[n];
@ -2608,6 +2600,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
//........................................................................
//..............carry out relaxation process..............................
//..........Toelke, Fruediger et. al. 2006................................
//---------------- NO higher-order force -------------------------------//
if (C == 0.0) nx = ny = nz = 0.0;
m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2);
@ -2632,6 +2625,43 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//----------------With higher-order force ------------------------------//
//if (C == 0.0) nx = ny = nz = 0.0;
//m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1)
// + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
//m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2)
// + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity;
//jx = jx + Fx;
//m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
//jy = jy + Fy;
//m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
//jz = jz + Fz;
//m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
//m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9)
// + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
////m10 = m10 + rlx_setA*( - m10);
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10)
// + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity;
//m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11)
// + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
////m12 = m12 + rlx_setA*( - m12);
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12)
// + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity;
//m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13);
// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
//m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14);
// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
//m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15);
// + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity;
//m16 = m16 + rlx_setB*( - m16);
//m17 = m17 + rlx_setB*( - m17);
//m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//.................inverse transformation......................................................
// q=0
@ -2748,6 +2778,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@ -2764,6 +2798,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
// Cq = {0,1,0}
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@ -2779,6 +2817,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
// Cq = {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@ -2790,6 +2832,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
Aq[6*Np+n] = a2;
Bq[6*Np+n] = b2;
//...............................................
}
}

43
cuda/GreyscaleColor.cu
View File

@ -1450,7 +1450,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
//CP: capillary penalty
// also turn off recoloring for grey nodes
__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, 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){
@ -1478,6 +1478,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
double Fx,Fy,Fz;
double Fcpx,Fcpy,Fcpz;//capillary penalty force
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
const double mrt_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802;
@ -1504,6 +1505,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
porosity = Poros[n];
perm = Perm[n];
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
@ -2165,6 +2168,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@ -2184,6 +2191,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
// Cq = {0,1,0}
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@ -2204,6 +2215,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
// Cq = {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@ -2226,7 +2241,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
//CP: capillary penalty
// also turn off recoloring for grey nodes
__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, 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){
int ijk,nn,n;
@ -2249,6 +2264,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
double Fx,Fy,Fz;
double Fcpx,Fcpy,Fcpz;//capillary penalty force
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
const double mrt_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802;
@ -2276,6 +2292,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
porosity = Poros[n];
perm = Perm[n];
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
@ -2870,6 +2888,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@ -2886,6 +2908,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
// Cq = {0,1,0}
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@ -2901,6 +2927,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
// Cq = {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@ -4500,12 +4530,13 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl
//Model-1 & 4 with capillary pressure penalty
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, Poros, Perm, Vel, Pressure,
dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm, Vel, Pressure,
rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AAeven_GreyscaleColor_CP: %s \n",cudaGetErrorString(err));
@ -4515,11 +4546,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
//Model-1 & 4 with capillary pressure penalty
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 *Poros,double *Perm,double *Vel,double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, Poros, Perm,Vel,Pressure,
dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm,Vel,Pressure,
rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
cudaError_t err = cudaGetLastError();

20
docs/Makefile Normal file
View File

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

19
docs/README.md Normal file
View File

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

@ -0,0 +1,164 @@
============
Installation
============
Sample scripts to configure and build LBPM are included in the `sample_scripts` directory. To build this package, MPI and cmake are required, along with a compiler that supports the C++-14 standard. Building in source is not supported.
The essential dependencies needed to build LBPM are:
1. cmake (version 3.9 or higher)
2. MPI
3. HDF5
4. silo
*************************
Building Dependencies
*************************
(skip if they are already installed on your system)
1. Download third-party library source files
.. code-block:: bash
zlib-1.2.11.tar.gz
hdf5-1.8.12.tar.gz
silo-4.10.2.tar.gz
2. Set up path where you want to install
.. code-block:: bash
export MPI_DIR=/path/to/mpi/
export LBPM_ZLIB_DIR=/path/to/zlib
export LBPM_HDF5_DIR=/path/to/hdf5
export LBPM_SILO_DIR=/path/to/silo
3. Build third-party library dependencies
.. code-block:: bash
tar -xzvf zlib-1.2.11.tar.gz
tar -xzvf hdf5-1.8.12.tar.gz
tar -xzvf silo-4.10.2.tar.gz
cd zlib-1.2.11
./configure --prefix=$LBPM_ZLIB_DIR && make && make install
cd ../hdf5-1.8.12
CC=$MPI_DIR/bin/mpicc CXX=$MPI_DIR/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \
./configure --prefix=$LBPM_HDF5_DIR --enable-parallel --enable-shared --with-zlib=$LBPM_ZLIB_DIR
make && make install
cd ../silo-4.10.2
CC=$MPI_DIR/bin/mpicc CXX=$MPI_DIR/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \
./configure --prefix=$LBPM_SILO_DIR -with-hdf5=$LBPM_HDF5_DIR/include,$LBPM_HDF5_DIR/lib --enable-static
make && make install
*************************
Building LBPM
*************************
Many HPC systems will include all of these dependencies, and LBPM can be built simply by setting the paths to the associated libraries in the cmake configure line.
The steps typically used to build LBPM are as follows:
1. Set environment variables for the source directory `$LBPM_WIA_SOURCE` and build directory `$LBPM_WIA_DIR`
.. code-block:: bash
export LBPM_SOURCE=/path/to/source/LBPM
export LBPM_DIR=/path/to/build/LBPM
2. Set environment variables for the path to HDF5 and SILO (required), and optionally for TimerUtility and NetCDF (optional)
.. code-block:: bash
export LBPM_HDF5_DIR=/path/to/hdf5
export LBPM_SILO_DIR=/path/to/silo
export LBPM_TIMER_DIR=/path/to/timer
export LBPM_NETCDF_DIR=/path/to/netcdf
3. Create the build directory and navigate to it
.. code-block:: bash
mkdir $LBPM_WIA_DIR
cd $LBPM_WIA_DIR
4. Configure the project. Numerous scripts exist to build LBPM on different HPC clusters, which are available in the `$LBPM_SOURCE/sample_scripts/` directory. It is often possible to run these scripts directly if one exists for a system similar to the one you are building on. For a standard CPU build:
.. code-block:: bash
cmake \
-D CMAKE_BUILD_TYPE:STRING=Release \
-D CMAKE_C_COMPILER:PATH=mpicc \
-D CMAKE_CXX_COMPILER:PATH=mpicxx \
-D CMAKE_C_FLAGS="-fPIC" \
-D CMAKE_CXX_FLAGS="-fPIC" \
-D CMAKE_CXX_STD=14 \
-D USE_TIMER=0 \
-D TIMER_DIRECTORY=$LBPM_TIMER_DIR \
-D USE_NETCDF=0 \
-D NETCDF_DIRECTORY=$LBPM_NETCDF_DIR \
-D USE_SILO=1 \
-D HDF5_DIRECTORY=$LBPM_HDF5_DIR \
-D SILO_DIRECTORY=$LBPM_SILO_DIR \
-D USE_CUDA=0 \
$LBPM_SOURCE
For GPU support, it is necessary to have CUDA along with a GPU-aware MPI implementation. Otherwise, the LBPM routines should behave identically irrespective of the underlying hardware.
.. code-block:: bash
cmake \
-D CMAKE_BUILD_TYPE:STRING=Release \
-D CMAKE_C_COMPILER:PATH=mpicc \
-D CMAKE_CXX_COMPILER:PATH=mpicxx \
-D CMAKE_C_FLAGS="-fPIC" \
-D CMAKE_CXX_FLAGS="-fPIC" \
-D CMAKE_CXX_STD=14 \
-D USE_TIMER=0 \
-D TIMER_DIRECTORY=$LBPM_TIMER_DIR \
-D USE_NETCDF=0 \
-D NETCDF_DIRECTORY=$LBPM_NETCDF_DIR \
-D USE_SILO=1 \
-D HDF5_DIRECTORY=$LBPM_HDF5_DIR \
-D SILO_DIRECTORY=$LBPM_SILO_DIR \
-D USE_CUDA=1 \
-D CMAKE_CUDA_FLAGS="-arch sm_70" \
$LBPM_SOURCE
5. Build the project (using four cores to build)
.. code-block:: bash
make -j4
6. Install the project
.. code-block:: bash
make install
7. Run the tests to make sure they execute correctly (on a cluster, it is recommended to run these using the batch system rather than on the head node)
.. code-block:: bash
ctest
*************************
Sample Scripts
*************************
The LBPM repository contains sample scripts showing successful CMake configuration, build and
install steps for a range of systems. Refer to the project sub-directory below for these examples.
.. code-block:: bash
ls $LBPM_SOURCE/sample_scripts

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

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

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

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

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

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

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

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

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

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

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

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

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

@ -0,0 +1,43 @@
======================================
Color model -- Centrifuge Protocol
======================================
The centrifuge protocol is designed to mimic SCAL centrifuge experiments that
are used to infer the capillary pressure. The LBPM centrifuge protocol is
constructed as an unsteady simulation with constant pressure boundary conditions
and zero pressure drop across the sample. This will enforce the following key values
* ``BC = 3`` -- constant pressure boundary condition
* ``din = 1.0`` -- inlet pressure value
* ``dout = 1.0`` -- outlet pressure value
By default LBPM will populate the inlet reservoir with fluid A (usually the non-wetting fluid)
and the outlet reservoir with fluid B (usually water). Flow is induced by setting an external
body force to generate displacement in the ``z`` direction. If the body force is set to
zero, e.g. ``F = 0, 0, 0``, the simulation will produce spontaneous imbibition, with the
balancing point being determined based on zero pressure drop across the sample. Setting
an external body force will shift the capillary pressure. Setting a positive force will
cause fluid A to be forced into the sample. Once steady conditions are achieved,
the pressure of fluid A will be larger than fluid B. Alternatively, if the driving force is
negative then fluid B will be forced into the sample, and the steady-state configuration
will stabilize to a configuration where fluid B has a larger pressure compared to fluid A.
The capillary pressure is thereby inferred based on the body force.
In a conventional SCAL experiment the centrifugal forces are proportional to the density
difference between fluids. While this is still true for LBPM simulation, the body force will
still be effective even if there is no difference in density between the fluids.
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).
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
.. code-block:: bash
protocol = "centrifuge"
F = 0, 0, -1.0e-5

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

177
example/DiscPack/DiscPack.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

247
hip/GreyscaleColor.cu
View File

@ -1450,9 +1450,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
//CP: capillary penalty
// also turn off recoloring for grey nodes
__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, 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, double W, int strideY, int strideZ, int start, int finish, int Np){
double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
int n,nn,ijk,nread;
int nr1,nr2,nr3,nr4,nr5,nr6;
@ -1462,8 +1462,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
double fq;
// conserved momemnts
double rho,jx,jy,jz;
//double vx,vy,vz,v_mag;
//double ux,uy,uz,u_mag;
double ux,uy,uz;
// non-conserved moments
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
@ -1473,18 +1471,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
double C,nx,ny,nz; //color gradient magnitude and direction
double phi,tau,rho0,rlx_setA,rlx_setB;
//double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002)
double porosity;
double perm;//voxel permeability
//double c0, c1; //Guo's model parameters
double tau_eff;
double mu_eff;//kinematic viscosity
double nx_gs,ny_gs,nz_gs;//grey-solid color gradient
double nx_phase,ny_phase,nz_phase,C_phase;
double Fx,Fy,Fz;
double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18;
double gp3,gp5,gp7;
double Fcpx,Fcpy,Fcpz;//capillary penalty force
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
const double mrt_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802;
@ -1510,9 +1504,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
porosity = Poros[n];
perm = Perm[n];
nx_gs = GreySolidGrad[n+0*Np];
ny_gs = GreySolidGrad[n+1*Np];
nz_gs = GreySolidGrad[n+2*Np];
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
@ -1534,98 +1528,61 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
gp1 = Psi[nn];
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
gp2 = Psi[nn];
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Phi[nn]; // get neighbor for phi - 3
gp3 = Psi[nn];
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Phi[nn]; // get neighbor for phi - 4
gp4 = Psi[nn];
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Phi[nn]; // get neighbor for phi - 5
gp5 = Psi[nn];
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Phi[nn]; // get neighbor for phi - 6
gp6 = Psi[nn];
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Phi[nn]; // get neighbor for phi - 7
gp7 = Psi[nn];
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Phi[nn]; // get neighbor for phi - 8
gp8 = Psi[nn];
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Phi[nn]; // get neighbor for phi - 9
gp9 = Psi[nn];
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Phi[nn]; // get neighbor for phi - 10
gp10 = Psi[nn];
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Phi[nn]; // get neighbor for phi - 11
gp11 = Psi[nn];
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Phi[nn]; // get neighbor for phi - 12
gp12 = Psi[nn];
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Phi[nn]; // get neighbor for phi - 13
gp13 = Psi[nn];
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Phi[nn]; // get neighbor for phi - 14
gp14 = Psi[nn];
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Phi[nn]; // get neighbor for phi - 15
gp15 = Psi[nn];
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Phi[nn]; // get neighbor for phi - 16
gp16 = Psi[nn];
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Phi[nn]; // get neighbor for phi - 17
gp17 = Psi[nn];
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Phi[nn]; // get neighbor for phi - 18
gp18 = Psi[nn];
//............Compute the Color Gradient...................................
nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase);
//correct the normal color gradient by considering the effect of grey solid
nx = nx_phase + (1.0-porosity)*nx_gs;
ny = ny_phase + (1.0-porosity)*ny_gs;
nz = nz_phase + (1.0-porosity)*nz_gs;
if (C_phase==0.0){
nx = nx_phase;
ny = ny_phase;
nz = nz_phase;
}
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//............Compute the Greyscale Potential Gradient.....................
// Fcpx = 0.0;
@ -1648,14 +1605,24 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
// //ny = Fcpy/Fcp_mag;
// //nz = Fcpz/Fcp_mag;
// }
Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
Fcpx = nx;
Fcpy = ny;
Fcpz = nz;
double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
if (Fcp_mag==0.0); Fcpx=Fcpy=Fcpz=0.0;
//NOTE for open node (porosity=1.0),Fcp=0.0
Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm);
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=0
fq = dist[n];
rho = fq;
@ -2201,6 +2168,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@ -2220,6 +2191,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
// Cq = {0,1,0}
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@ -2240,6 +2215,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
// Cq = {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@ -2262,15 +2241,13 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
//CP: capillary penalty
// also turn off recoloring for grey nodes
__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, 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, double W, int strideY, int strideZ, int start, int finish, int Np){
double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
int ijk,nn,n;
double fq;
// conserved momemnts
double rho,jx,jy,jz;
//double vx,vy,vz,v_mag;
//double ux,uy,uz,u_mag;
double ux,uy,uz;
// non-conserved moments
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
@ -2280,18 +2257,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
double C,nx,ny,nz; //color gradient magnitude and direction
double phi,tau,rho0,rlx_setA,rlx_setB;
//double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002)
double porosity;
double perm;//voxel permeability
//double c0, c1; //Guo's model parameters
double tau_eff;
double mu_eff;//kinematic viscosity
double nx_gs,ny_gs,nz_gs;//grey-solid color gradient
double nx_phase,ny_phase,nz_phase,C_phase;
double Fx,Fy,Fz;
double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18;
double gp3,gp5,gp7;
double Fcpx,Fcpy,Fcpz;//capillary penalty force
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
const double mrt_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802;
@ -2315,11 +2288,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
// read the component number densities
nA = Den[n];
nB = Den[Np + n];
porosity = Poros[n];
perm = Perm[n];
nx_gs = GreySolidGrad[n+0*Np];
ny_gs = GreySolidGrad[n+1*Np];
nz_gs = GreySolidGrad[n+2*Np];
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
@ -2341,98 +2315,61 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
gp1 = Psi[nn];
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
gp2 = Psi[nn];
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Phi[nn]; // get neighbor for phi - 3
gp3 = Psi[nn];
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Phi[nn]; // get neighbor for phi - 4
gp4 = Psi[nn];
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Phi[nn]; // get neighbor for phi - 5
gp5 = Psi[nn];
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Phi[nn]; // get neighbor for phi - 6
gp6 = Psi[nn];
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Phi[nn]; // get neighbor for phi - 7
gp7 = Psi[nn];
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Phi[nn]; // get neighbor for phi - 8
gp8 = Psi[nn];
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Phi[nn]; // get neighbor for phi - 9
gp9 = Psi[nn];
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Phi[nn]; // get neighbor for phi - 10
gp10 = Psi[nn];
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Phi[nn]; // get neighbor for phi - 11
gp11 = Psi[nn];
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Phi[nn]; // get neighbor for phi - 12
gp12 = Psi[nn];
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Phi[nn]; // get neighbor for phi - 13
gp13 = Psi[nn];
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Phi[nn]; // get neighbor for phi - 14
gp14 = Psi[nn];
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Phi[nn]; // get neighbor for phi - 15
gp15 = Psi[nn];
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Phi[nn]; // get neighbor for phi - 16
gp16 = Psi[nn];
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Phi[nn]; // get neighbor for phi - 17
gp17 = Psi[nn];
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Phi[nn]; // get neighbor for phi - 18
gp18 = Psi[nn];
//............Compute the Color Gradient...................................
nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase);
//correct the normal color gradient by considering the effect of grey solid
nx = nx_phase + (1.0-porosity)*nx_gs;
ny = ny_phase + (1.0-porosity)*ny_gs;
nz = nz_phase + (1.0-porosity)*nz_gs;
if (C_phase==0.0){
nx = nx_phase;
ny = ny_phase;
nz = nz_phase;
}
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//............Compute the Greyscale Potential Gradient.....................
// Fcpx = 0.0;
@ -2455,14 +2392,23 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
// ny = Fcpy/Fcp_mag;
// nz = Fcpz/Fcp_mag;
// }
Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
Fcpx = nx;
Fcpy = ny;
Fcpz = nz;
double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
if (Fcp_mag==0.0); Fcpx=Fcpy=Fcpz=0.0;
//NOTE for open node (porosity=1.0),Fcp=0.0
Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm);
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=0
fq = dist[n];
@ -2942,6 +2888,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
@ -2958,6 +2908,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
// Cq = {0,1,0}
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
@ -2973,6 +2927,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
// Cq = {0,0,1}
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
@ -2989,6 +2947,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
}
}
__global__ void dvc_ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np){
int idx;
double nA,nB;
@ -4572,12 +4531,12 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl
//Model-1 & 4 with capillary pressure penalty
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi,double *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, Psi, GreySolidGrad, Poros, Perm, Vel, Pressure,
rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, W, strideY, strideZ, start, finish, Np);
dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm, Vel, Pressure,
rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("hip error in ScaLBL_D3Q19_AAeven_GreyscaleColor_CP: %s \n",hipGetErrorString(err));
@ -4587,12 +4546,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
//Model-1 & 4 with capillary pressure penalty
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, Psi, GreySolidGrad, Poros, Perm,Vel,Pressure,
rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, W, strideY, strideZ, start, finish, Np);
dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm,Vel,Pressure,
rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
@ -4600,51 +4559,3 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *M
}
}
//extern "C" void ScaLBL_Update_GreyscalePotential(int *Map, double *Phi, double *Psi, double *Poro, double *Perm, double alpha, double W,
// int start, int finish, int Np){
//
// dvc_ScaLBL_Update_GreyscalePotential<<<NBLOCKS,NTHREADS >>>(Map, Phi, Psi, Poro, Perm, alpha, W, start, finish, Np);
//
// hipError_t err = hipGetLastError();
// if (hipSuccess != err){
// printf("hip error in ScaLBL_Update_GreyscalePotential: %s \n",hipGetErrorString(err));
// }
//}
////Model-2&3
//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den,
// double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel,
// double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
//
// //cudaProfilerStart();
// //cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor, cudaFuncCachePreferL1);
//
// dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm, Vel,
// rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np);
// hipError_t err = hipGetLastError();
// if (hipSuccess != err){
// printf("hip error in ScaLBL_D3Q19_AAeven_GreyscaleColor: %s \n",hipGetErrorString(err));
// }
// //cudaProfilerStop();
//
//}
//
////Model-2&3
//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
// double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,
// double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
//
// //cudaProfilerStart();
// //cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor, cudaFuncCachePreferL1);
//
// dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm,Vel,
// rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np);
//
// hipError_t err = hipGetLastError();
// if (hipSuccess != err){
// printf("hip error in ScaLBL_D3Q19_AAodd_GreyscaleColor: %s \n",hipGetErrorString(err));
// }
// //cudaProfilerStop();
//}

947
models/ColorModel.cpp
View File

File diff suppressed because it is too large Load Diff

8
models/ColorModel.h
View File

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

73
models/FreeLeeModel.cpp
View File

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

33
models/GreyscaleColorModel.cpp
View File

@ -309,18 +309,26 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
// oil-wet < 0
// neutral = 0 (i.e. no penalty)
double *GreySolidW_host = new double [Np];
double *GreySn_host = new double [Np];
double *GreySw_host = new double [Np];
size_t NLABELS=0;
signed char VALUE=0;
double AFFINITY=0.f;
double Sn,Sw;//end-point saturation of greynodes set by users
auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
auto SnList = greyscaleColor_db->getVector<double>( "GreySnList" );
auto SwList = greyscaleColor_db->getVector<double>( "GreySwList" );
NLABELS=LabelList.size();
if (NLABELS != AffinityList.size()){
ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
}
if (NLABELS != SnList.size() || NLABELS != SwList.size()){
ERROR("Error: GreySolidLabels, GreySnList, and GreySwList must be the same length! \n");
}
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
@ -328,16 +336,22 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
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
Sn=99.0;
Sw=-99.0;
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
Sn = SnList[idx];
Sw = SwList[idx];
idx = NLABELS;
}
}
int idx = Map(i,j,k);
if (!(idx < 0)){
GreySolidW_host[idx] = AFFINITY;
GreySn_host[idx] = Sn;
GreySw_host[idx] = Sw;
}
}
}
@ -348,15 +362,22 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
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);
Sn=SnList[idx];
Sw=SwList[idx];
//printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
printf(" grey-solid label=%d, grey-solid affinity=%.3g, grey-solid Sn=%.3g, grey-solid Sw=%.3g\n",VALUE,AFFINITY,Sn,Sw);
}
printf("NOTE: grey-solid affinity>0: water-wet || grey-solid affinity<0: oil-wet \n");
}
ScaLBL_CopyToDevice(GreySolidW, GreySolidW_host, Np*sizeof(double));
ScaLBL_CopyToDevice(GreySn, GreySn_host, Np*sizeof(double));
ScaLBL_CopyToDevice(GreySw, GreySw_host, Np*sizeof(double));
ScaLBL_Comm->Barrier();
delete [] GreySolidW_host;
delete [] GreySn_host;
delete [] GreySw_host;
}
////----------------------------------------------------------------------------------------------------------//
@ -598,6 +619,8 @@ void ScaLBL_GreyscaleColorModel::Create(){
//ScaLBL_AllocateDeviceMemory((void **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz);
//ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreySolidW, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreySn, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreySw, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Porosity_dvc, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Permeability_dvc, sizeof(double)*Np);
//...........................................................................
@ -921,7 +944,7 @@ 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,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,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
@ -950,7 +973,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,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
@ -983,7 +1006,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
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,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,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
@ -1012,7 +1035,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,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

2
models/GreyscaleColorModel.h
View File

@ -68,6 +68,8 @@ public:
//double *GreySolidPhi; //Model 2 & 3
//double *GreySolidGrad;//Model 1 & 4
double *GreySolidW;
double *GreySn;
double *GreySw;
//double *ColorGrad;
double *Velocity;
double *Pressure;

49
models/GreyscaleModel.cpp
View File

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

74
models/IonModel.cpp
View File

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

4
models/IonModel.h
View File

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

8
models/StokesModel.cpp
View File

@ -59,6 +59,7 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
BoundaryCondition = 0;
if (stokes_db->keyExists( "BC" )){
BoundaryCondition = stokes_db->getScalar<int>( "BC" );
}
if (stokes_db->keyExists( "tolerance" )){
tolerance = stokes_db->getScalar<double>( "tolerance" );
@ -110,6 +111,7 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
void ScaLBL_StokesModel::ReadParams(string filename){
//NOTE the max time step is left unspecified
// read the input database
db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" );
@ -298,8 +300,10 @@ void ScaLBL_StokesModel::AssignZetaPotentialSolid(double *zeta_potential_solid)
ERROR("Error: LB Single-Fluid Solver: SolidLabels and ZetaPotentialSolidList must be the same length! \n");
}
double label_count[NLABELS];
double label_count_global[NLABELS];
double *label_count;
double *label_count_global;
label_count = new double [NLABELS];
label_count_global = new double [NLABELS];
for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;

67
tests/DataAggregator.cpp
View File

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

14
tests/GenerateSphereTest.cpp
View File

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

225
tests/lbpm_color_simulator.cpp
View File

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

1
tests/lbpm_freelee_SingleFluidBGK_simulator.cpp
View File

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

1
tests/lbpm_freelee_simulator.cpp
View File

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