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

Merge branch 'master' into FOM

Browse Source
This commit is contained in:
James E McClure 2021-06-16 16:24:13 -04:00
commit 399191f8dc
75 changed files with 15972 additions and 777 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists.txt
View File

@ -1,6 +1,8 @@
# Set some CMake properties
CMAKE_MINIMUM_REQUIRED( VERSION 3.9 )
CMAKE_POLICY( SET CMP0115 OLD )
if( ${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.20.0")
CMAKE_POLICY( SET CMP0115 OLD )
endif()
MESSAGE("====================")
@ -16,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}" )

2
IO/IOHelpers.h
View File

@ -57,6 +57,6 @@ inline std::vector<std::string> splitList( const char *line, const char token )
}
}; // namespace IO
} // namespace IO
#endif

7
IO/silo.h
View File

@ -18,7 +18,8 @@ typedef int DBfile;
#endif
namespace IO::silo {
namespace IO {
namespace silo {
enum FileMode { READ, WRITE, CREATE };
@ -256,7 +257,9 @@ void writeMultiVar( DBfile *fid, const std::string &varname,
const std::vector<std::string> &subVarNames, const std::vector<int> &subVarTypes );
}; // namespace IO::silo
} // namespace silo
} // namespace IO
#endif
#include "IO/silo.hpp"

6
IO/silo.hpp
View File

@ -13,7 +13,8 @@
#include <silo.h>
namespace IO::silo {
namespace IO {
namespace silo {
/****************************************************
@ -413,7 +414,8 @@ Array<TYPE> readTriMeshVariable( DBfile *fid, const std::string &varname )
}
}; // namespace IO::silo
} // namespace silo
} // namespace IO
#endif

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);

89
analysis/SubPhase.cpp
View File

@ -93,7 +93,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 pw pn wet\n");
fprintf(TIMELOG,"sw krw krn krwf krnf vw vn force pw pn wet\n");
}
}
}
@ -167,7 +167,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 +245,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 +288,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 +300,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 +335,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 +357,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 +367,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\n",saturation,krw,krn,h*water_flow_rate,h*not_water_flow_rate, 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 +397,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 +436,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 +447,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);
}
@ -606,8 +655,8 @@ void SubPhase::Full(){
double uy = Vel_y(n);
double uz = Vel_z(n);
if (DelPhi(n) > 1e-3){
// interface region
if (DelPhi(n) > 1e-3 && SDs(n) < 3.0 ){
// film region
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));

47
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)
@ -320,6 +317,8 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
DoubleArray phase(nx,ny,nz);
IntArray phase_label(nx,ny,nz);
Array<char> ID(nx,ny,nz);
fillHalo<char> fillChar(Dm->Comm,Dm->rank_info,{nx-2,ny-2,nz-2},{1,1,1},0,1);
int n;
double final_void_fraction;
@ -337,10 +336,11 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
count += 1.0;
id[n] = 2;
}
ID(i,j,k) = id[n];
}
}
}
fillChar.fill(ID);
Dm->Comm.barrier();
// total Global is the number of nodes in the pore-space
@ -351,7 +351,8 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
if (rank==0) printf("Volume fraction for morphological opening: %f \n",volume_fraction);
if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);
// Communication buffers
/* // Communication buffers
signed char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
signed char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
signed char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
@ -400,8 +401,9 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
//......................................................................................
int sendtag,recvtag;
sendtag = recvtag = 7;
*/
int ii,jj,kk;
int imin,jmin,kmin,imax,jmax,kmax;
int Nx = nx;
int Ny = ny;
int Nz = nz;
@ -413,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);
@ -453,11 +452,11 @@ 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[nn] == 2 && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
if (ID(ii,jj,kk) == 2 && dsq <= (Rcrit_new+1)*(Rcrit_new+1)){
LocalNumber+=1.0;
id[nn]=1;
//id[nn]=1;
ID(ii,jj,kk)=1;
}
}
}
@ -468,8 +467,9 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
}
}
}
fillChar.fill(ID);
// Pack and send the updated ID values
PackID(Dm->sendList("x"), Dm->sendCount("x") ,sendID_x, id);
/* PackID(Dm->sendList("x"), Dm->sendCount("x") ,sendID_x, id);
PackID(Dm->sendList("X"), Dm->sendCount("X") ,sendID_X, id);
PackID(Dm->sendList("y"), Dm->sendCount("y") ,sendID_y, id);
PackID(Dm->sendList("Y"), Dm->sendCount("Y") ,sendID_Y, id);
@ -527,12 +527,12 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
UnpackID(Dm->recvList("YZ"), Dm->recvCount("YZ") ,recvID_YZ, id);
//......................................................................................
// double GlobalNumber = Dm->Comm.sumReduce( LocalNumber );
*/
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
if (id[n] == 1){
if (ID(i,j,k) == 1){
phase(i,j,k) = 1.0;
}
else
@ -550,9 +550,10 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
if (id[n] == 1 && phase_label(i,j,k) > 1){
id[n] = 2;
if (ID(i,j,k) == 1 && phase_label(i,j,k) > 1){
ID(i,j,k) = 2;
}
id[n] = ID(i,j,k);
}
}
}
@ -571,7 +572,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
// nwp
phase(i,j,k) = -1.0;
}
else{
else{i
// treat solid as WP since films can connect
phase(i,j,k) = 1.0;
}

1
analysis/runAnalysis.cpp
View File

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

6
cmake/macros.cmake
View File

@ -1,6 +1,8 @@
CMAKE_MINIMUM_REQUIRED(VERSION 3.18.3)
CMAKE_MINIMUM_REQUIRED(VERSION 3.10.0)
CMAKE_POLICY( SET CMP0057 NEW )
CMAKE_POLICY( SET CMP0115 OLD )
if( ${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.20.0")
CMAKE_POLICY( SET CMP0115 OLD )
endif()
INCLUDE(CheckCCompilerFlag)
INCLUDE(CheckCSourceCompiles)

17
common/Domain.cpp
View File

@ -138,7 +138,7 @@ void Domain::initialize( std::shared_ptr<Database> db )
if (rank_info.kz < nproc[2]-1) outlet_layers_z = 0;
// Fill remaining variables
N = Nx*Ny*Nz;
Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0;
Volume = nx*ny*nz*nproc[0]*nproc[1]*nproc[2]*1.0;
if (myrank==0) printf("voxel length = %f micron \n", voxel_length);
@ -620,12 +620,16 @@ void Domain::Decomp( const std::string& Filename )
Comm.recv(id.data(),N,0,15);
}
Comm.barrier();
ComputePorosity();
delete [] SegData;
}
void Domain::ComputePorosity(){
// Compute the porosity
double sum;
double sum_local=0.0;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocx()*nprocy()*nprocz());
if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
//.........................................................
for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-1;k++){
for (int j=1;j<Ny-1;j++){
@ -640,8 +644,8 @@ void Domain::Decomp( const std::string& Filename )
sum = Comm.sumReduce(sum_local);
porosity = sum*iVol_global;
if (rank()==0) printf("Media porosity = %f \n",porosity);
//.........................................................
delete [] SegData;
//.........................................................
}
void Domain::AggregateLabels( const std::string& filename ){
@ -1450,4 +1454,3 @@ void Domain::AggregateLabels( const std::string& filename, DoubleArray &UserData
Comm.barrier();
}

1
common/Domain.h
View File

@ -166,6 +166,7 @@ public: // Public variables (need to create accessors instead)
std::vector<signed char> id;
void ReadIDs();
void ComputePorosity();
void Decomp( const std::string& filename );
void CommunicateMeshHalo(DoubleArray &Mesh);
void CommInit();

168
common/ScaLBL.cpp
View File

@ -973,7 +973,7 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis
}
void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np)
void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC)
{
int idx,i,j,k;
@ -1046,12 +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;
/* 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++){
@ -1064,36 +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){
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){
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){
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){
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){
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){
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;
}
}
@ -1111,72 +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){
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){
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){
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){
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){
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){
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){
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){
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){
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){
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){
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){
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;
}
}
@ -1186,10 +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));
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;
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){
@ -1204,6 +1356,14 @@ void ScaLBL_Communicator::SolidNeumannD3Q7(double *fq, double *BoundaryValue){
ScaLBL_Solid_Neumann_D3Q7(fq,BoundaryValue,bb_dist,bb_interactions,n_bb_d3q7);
}
void ScaLBL_Communicator::SolidSlippingVelocityBCD3Q19(double *fq, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0, double den_scale,double h, double time_conv){
// fq is a D3Q19 distribution
// BoundaryValues is a list of values to assign at bounce-back solid sites
ScaLBL_Solid_SlippingVelocityBC_D3Q19(fq,zeta_potential,ElectricField,SolidGrad,epsilon_LB,tau,rho0,den_scale,h,time_conv,
bb_dist,bb_interactions,fluid_boundary,lattice_weight,lattice_cx,lattice_cy,lattice_cz,n_bb_d3q19,N);
}
void ScaLBL_Communicator::SendD3Q19AA(double *dist){
// NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2

51
common/ScaLBL.h
View File

@ -87,6 +87,19 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *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 *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_Update_GreyscalePotential(int *Map, double *Phi, double *Psi, double *Poro, double *Perm, double alpha, double W,
// int start, int finish, int Np);
// ION TRANSPORT MODEL
extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *dist, double *Den, int start, int finish, int Np);
@ -193,11 +206,12 @@ extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double
extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad,
double rhonA, double rhoB, double tauM, double W, int start, int finish, int Np);
//extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi,
// double rhoA, double rhoB, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi,
double rhoA, double rhoB, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi,
double rhoA, double rhoB, int start, int finish, int Np);
//extern "C" void ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(int *Map, double *hq, double *Den, double *Phi,
// double rhoA, double rhoB, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(int *neighborList, int *Map, double *hq, double *Den, double *Phi, double *ColorGrad, double *Vel,
double rhoA, double rhoB, double tauM, double W, int start, int finish, int Np);
@ -214,6 +228,22 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double
double rhoA, double rhoB, double tauA, double tauB, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined(int *neighborList, int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined_HigherOrder(int *neighborList, int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined_HigherOrder(int *Map, double *dist, double *hq, double *Den, double *Phi, double *mu_phi, double *Vel, double *Pressure, double *ColorGrad,
double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz,
int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure,
double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np);
@ -258,6 +288,12 @@ extern "C" void ScaLBL_Solid_Dirichlet_D3Q7(double *dist,double *BoundaryValue,i
extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int *BounceBackDist_list,int *BounceBackSolid_list,int N);
extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
int count, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np);
@ -337,9 +373,11 @@ public:
void RecvHalo(double *data);
void RecvGrad(double *Phi, double *Gradient);
void RegularLayout(IntArray map, const double *data, DoubleArray &regdata);
void SetupBounceBackList(IntArray &Map, signed char *id, int Np);
void SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC=false);
void SolidDirichletD3Q7(double *fq, double *BoundaryValue);
void SolidNeumannD3Q7(double *fq, double *BoundaryValue);
void SolidSlippingVelocityBCD3Q19(double *fq, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epslion_LB, double tau, double rho0, double den_scale,double h, double time_conv);
// Routines to set boundary conditions
void Color_BC_z(int *Map, double *Phi, double *Den, double vA, double vB);
@ -414,6 +452,9 @@ private:
//......................................................................................
int *bb_dist;
int *bb_interactions;
int *fluid_boundary;
double *lattice_weight;
float *lattice_cx, *lattice_cy, *lattice_cz;
//......................................................................................
};

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;

51
cpu/D3Q7BC.cpp
View File

@ -30,6 +30,57 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int
}
}
extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
int count, int Np){
int idx;
int iq,ib,ifluidBC;
double value_b,value_q;
double Ex,Ey,Ez;
double Etx,Ety,Etz;//tangential part of electric field
double E_mag_normal;
double nsx,nsy,nsz;//unit normal solid gradient
double ubx,uby,ubz;//slipping velocity at fluid boundary nodes
float cx,cy,cz;//lattice velocity (D3Q19)
double LB_weight;//lattice weighting coefficient (D3Q19)
double cs2_inv = 3.0;//inverse of cs^2 for D3Q19
double nu_LB = (tau-0.5)/cs2_inv;
for (idx=0; idx<count; idx++){
iq = BounceBackDist_list[idx];
ib = BounceBackSolid_list[idx];
ifluidBC = FluidBoundary_list[idx];
value_b = zeta_potential[ib];//get zeta potential from a solid site
value_q = dist[iq];
//Load electric field and compute its tangential componet
Ex = ElectricField[ifluidBC+0*Np];
Ey = ElectricField[ifluidBC+1*Np];
Ez = ElectricField[ifluidBC+2*Np];
nsx = SolidGrad[ifluidBC+0*Np];
nsy = SolidGrad[ifluidBC+1*Np];
nsz = SolidGrad[ifluidBC+2*Np];
E_mag_normal = Ex*nsx+Ey*nsy+Ez*nsz;//magnitude of electric field in the direction normal to solid nodes
//compute tangential electric field
Etx = Ex - E_mag_normal*nsx;
Ety = Ey - E_mag_normal*nsy;
Etz = Ez - E_mag_normal*nsz;
ubx = -epsilon_LB*value_b*Etx/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
uby = -epsilon_LB*value_b*Ety/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
ubz = -epsilon_LB*value_b*Etz/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
//compute bounce-back distribution
LB_weight = lattice_weight[idx];
cx = lattice_cx[idx];
cy = lattice_cy[idx];
cz = lattice_cz[idx];
dist[iq] = value_q - 2.0*LB_weight*rho0*cs2_inv*(cx*ubx+cy*uby+cz*ubz);
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){
for (int idx=0; idx<count; idx++){
int n = list[idx];

2703
cpu/FreeLee.cpp
View File

File diff suppressed because it is too large Load Diff

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

1499
cpu/GreyscaleColor.cpp
View File

File diff suppressed because it is too large Load Diff

2
cpu/MixedGradient.cpp
View File

@ -7,7 +7,7 @@ extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradie
{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1},
{0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}};
int i,j,k,n,N;
int i,j,k,n;
int np,np2,nm; // neighbors
double v,vp,vp2,vm; // values at neighbors
double grad;

68
cuda/D3Q7BC.cu
View File

@ -38,6 +38,57 @@ __global__ void dvc_ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValu
}
}
__global__ void dvc_ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
int count, int Np)
{
int idx;
int iq,ib,ifluidBC;
double value_b,value_q;
double Ex,Ey,Ez;
double Etx,Ety,Etz;//tangential part of electric field
double E_mag_normal;
double nsx,nsy,nsz;//unit normal solid gradient
double ubx,uby,ubz;//slipping velocity at fluid boundary nodes
float cx,cy,cz;//lattice velocity (D3Q19)
double LB_weight;//lattice weighting coefficient (D3Q19)
double cs2_inv = 3.0;//inverse of cs^2 for D3Q19
double nu_LB = (tau-0.5)/cs2_inv;
idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
iq = BounceBackDist_list[idx];
ib = BounceBackSolid_list[idx];
ifluidBC = FluidBoundary_list[idx];
value_b = zeta_potential[ib];//get zeta potential from a solid site
value_q = dist[iq];
//Load electric field and compute its tangential componet
Ex = ElectricField[ifluidBC+0*Np];
Ey = ElectricField[ifluidBC+1*Np];
Ez = ElectricField[ifluidBC+2*Np];
nsx = SolidGrad[ifluidBC+0*Np];
nsy = SolidGrad[ifluidBC+1*Np];
nsz = SolidGrad[ifluidBC+2*Np];
E_mag_normal = Ex*nsx+Ey*nsy+Ez*nsz;//magnitude of electric field in the direction normal to solid nodes
//compute tangential electric field
Etx = Ex - E_mag_normal*nsx;
Ety = Ey - E_mag_normal*nsy;
Etz = Ez - E_mag_normal*nsz;
ubx = -epsilon_LB*value_b*Etx/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
uby = -epsilon_LB*value_b*Ety/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
ubz = -epsilon_LB*value_b*Etz/(nu_LB*rho0)*time_conv*time_conv/(h*h*1.0e-12)/den_scale;
//compute bounce-back distribution
LB_weight = lattice_weight[idx];
cx = lattice_cx[idx];
cy = lattice_cy[idx];
cz = lattice_cz[idx];
dist[iq] = value_q - 2.0*LB_weight*rho0*cs2_inv*(cx*ubx+cy*uby+cz*ubz);
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np)
{
int idx,n;
@ -410,6 +461,23 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValue, i
}
}
extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad,
double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv,
int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list,
double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
int count, int Np){
int GRID = count / 512 + 1;
dvc_ScaLBL_Solid_SlippingVelocityBC_D3Q19<<<GRID,512>>>(dist, zeta_potential, ElectricField, SolidGrad,
epsilon_LB, tau, rho0, den_scale, h, time_conv,
BounceBackDist_list, BounceBackSolid_list, FluidBoundary_list,
lattice_weight, lattice_cx, lattice_cy, lattice_cz,
count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_Solid_SlippingVelocityBC_D3Q19 (kernel): %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z<<<GRID,512>>>(list, dist, Vin, count, Np);

4358
cuda/FreeLee.cu
View File

File diff suppressed because it is too large Load Diff

1573
cuda/GreyscaleColor.cu
View File

File diff suppressed because it is too large Load Diff

20
docs/Makefile Normal file
View File

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

19
docs/README.md Normal file
View File

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

16
example/Bubble/input.db
View File

@ -12,6 +12,22 @@ Color {
ComponentAffinity = -1.0, 1.0, -1.0
}
FreeLee {
tauA = 1.0;
tauB = 1.0;
tauM = 1.0;//relaxation parameter for the phase field
rhoA = 1.0;
rhoB = 1.0;
gamma = 1.0e-4;//surface tension parameter in Lee model
W = 3.0; //theoretical interfacial thickness in Lee model; unit:[voxel]
F = 0, 0, 0
Restart = false
timestepMax = 1000
flux = 0.0
ComponentLabels = 0
ComponentAffinity = -1.0
}
Domain {
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 80, 80, 80 // Size of local domain (Nx,Ny,Nz)

1468
hip/FreeLee.cu
View File

File diff suppressed because it is too large Load Diff

1601
hip/GreyscaleColor.cu
View File

File diff suppressed because it is too large Load Diff

924
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:

176
models/FreeLeeModel.cpp
View File

@ -59,6 +59,28 @@ void ScaLBL_FreeLeeModel::getVelocity(DoubleArray &Vel_x, DoubleArray &Vel_y, Do
ScaLBL_Comm->Barrier(); comm.barrier();
}
void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray &regdata){
// Gets data (in optimized layout) from the HOST and stores in regular layout
// Primarly for debugging
int i,j,k,idx;
// initialize the array
regdata.fill(0.f);
double value;
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
idx=Map(i,j,k);
if (!(idx<0)){
value=data[idx];
regdata(i,j,k)=value;
}
}
}
}
}
void ScaLBL_FreeLeeModel::ReadParams(string filename){
// read the input database
db = std::make_shared<Database>( filename );
@ -77,8 +99,10 @@ void ScaLBL_FreeLeeModel::ReadParams(string filename){
Fx = Fy = Fz = 0.0;
gamma=1e-3;//surface tension
W=5.0;//interfacial thickness
beta = 12.0*gamma/W;
kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
//beta = 12.0*gamma/W;
//kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
beta = 0.75*gamma/W;
kappa = 0.375*gamma*W;//beta and kappa are related to surface tension \gamma
Restart=false;
din=dout=1.0;
flux=0.0;
@ -136,8 +160,10 @@ void ScaLBL_FreeLeeModel::ReadParams(string filename){
outletA=0.f;
outletB=1.f;
//update secondary parameters
beta = 12.0*gamma/W;
kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
//beta = 12.0*gamma/W;
//kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
beta = 0.75*gamma/W;
kappa = 0.375*gamma*W;//beta and kappa are related to surface tension \gamma
//if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details)
BoundaryCondition = 0;
@ -386,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;
@ -570,27 +598,28 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels_ChemPotential_ColorGrad()
DoubleArray PhaseField(Nx,Ny,Nz);
FILE *OUTFILE;
ScaLBL_Comm->RegularLayout(Map,mu_phi_host,PhaseField);
getData_RegularLayout(mu_phi_host,PhaseField);
sprintf(LocalRankFilename,"Chem_Init.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE);
ScaLBL_Comm->RegularLayout(Map,&ColorGrad_host[0],PhaseField);
getData_RegularLayout(&ColorGrad_host[0],PhaseField);
FILE *CGX_FILE;
sprintf(LocalRankFilename,"Gradient_X_Init.%05i.raw",rank);
CGX_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CGX_FILE);
fclose(CGX_FILE);
ScaLBL_Comm->RegularLayout(Map,&ColorGrad_host[Np],PhaseField);
getData_RegularLayout(&ColorGrad_host[Np],PhaseField);
FILE *CGY_FILE;
sprintf(LocalRankFilename,"Gradient_Y_Init.%05i.raw",rank);
CGY_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CGY_FILE);
fclose(CGY_FILE);
ScaLBL_Comm->RegularLayout(Map,&ColorGrad_host[2*Np],PhaseField);
getData_RegularLayout(&ColorGrad_host[2*Np],PhaseField);
FILE *CGZ_FILE;
sprintf(LocalRankFilename,"Gradient_Z_Init.%05i.raw",rank);
CGZ_FILE = fopen(LocalRankFilename,"wb");
@ -709,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);
@ -795,27 +759,31 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
// Compute the Phase indicator field
// Read for hq happens in this routine (requires communication)
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(NeighborList, dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(NeighborList, dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastExterior(), Np);
//ScaLBL_D3Q7_AAodd_FreeLee_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
// Halo exchange for phase field
ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_Comm_WideHalo->Recv(Phi);
//ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_Comm_WideHalo->Send(Phi);
//ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL
if (BoundaryCondition > 0 && BoundaryCondition < 5){
//TODO to be revised
// Need to add BC for hq!!!
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// Halo exchange for phase field
ScaLBL_Comm_WideHalo->Send(Phi);
//ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined(NeighborList, dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
// Set BCs
@ -832,7 +800,9 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(gqbar);
}
ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
//ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q19_AAodd_FreeLeeModel_Combined(NeighborList, dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier();
@ -841,24 +811,29 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
timestep++;
// Compute the Phase indicator field
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMA
ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastExterior(), Np);
//ScaLBL_D3Q7_AAeven_FreeLee_PhaseField(dvcMap, hq, Den, Phi, ColorGrad, Velocity, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
// Halo exchange for phase field
ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_Comm_WideHalo->Recv(Phi);
//ScaLBL_D3Q7_ComputePhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_Comm_WideHalo->Send(Phi);
//ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL
if (BoundaryCondition > 0 && BoundaryCondition < 5){
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// Halo exchange for phase field
ScaLBL_Comm_WideHalo->Send(Phi);
//ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
// Set boundary conditions
@ -874,7 +849,9 @@ double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){
ScaLBL_Comm->D3Q19_Reflection_BC_z(gqbar);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(gqbar);
}
ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
//ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB,
// kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q19_AAeven_FreeLeeModel_Combined(dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, ColorGrad, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nxh, Nxh*Nyh, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier();
//************************************************************************
@ -1055,6 +1032,13 @@ void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){
fwrite(PhaseField.data(),8,N,PFILE);
fclose(PFILE);
ScaLBL_Comm->RegularLayout(Map,mu_phi,PhaseField);
FILE *CHEMFILE;
sprintf(LocalRankFilename,"ChemPotential.%05i.raw",rank);
CHEMFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,CHEMFILE);
fclose(CHEMFILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
FILE *VELX_FILE;
sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);

1
models/FreeLeeModel.h
View File

@ -84,6 +84,7 @@ public:
void getPhase(DoubleArray &PhaseValues);
void getPotential(DoubleArray &PressureValues, DoubleArray &MuValues);
void getVelocity(DoubleArray &Vx, DoubleArray &Vy, DoubleArray &Vz);
void getData_RegularLayout(const double *data, DoubleArray &regdata);
DoubleArray SignDist;

478
models/GreyscaleColorModel.cpp
View File

@ -17,7 +17,7 @@ void DeleteArray( const TYPE *p )
ScaLBL_GreyscaleColorModel::ScaLBL_GreyscaleColorModel(int RANK, int NP, const Utilities::MPI& COMM):
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),tauA_eff(0),tauB_eff(0),rhoA(0),rhoB(0),alpha(0),beta(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0),GreyPorosity(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0),GreyPorosity(0),RecoloringOff(0),
Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{
REVERSE_FLOW_DIRECTION = false;
@ -43,6 +43,8 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
Restart=false;
din=dout=1.0;
flux=0.0;
RecoloringOff = false;
//W=1.0;
// Color Model parameters
if (greyscaleColor_db->keyExists( "timestepMax" )){
@ -85,6 +87,9 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){
if (greyscaleColor_db->keyExists( "flux" )){
flux = greyscaleColor_db->getScalar<double>( "flux" );
}
if (greyscaleColor_db->keyExists( "RecoloringOff" )){
RecoloringOff = greyscaleColor_db->getScalar<bool>( "RecoloringOff" );
}
inletA=1.f;
inletB=0.f;
outletA=0.f;
@ -212,6 +217,7 @@ void ScaLBL_GreyscaleColorModel::ReadInput(){
}
void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
{
// Initialize impermeability solid nodes and grey nodes
@ -256,6 +262,7 @@ void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
label_count[idx] += 1.0;
idx = NLABELS;
//Mask->id[n] = 0; // set mask to zero since this is an immobile component
@ -290,31 +297,38 @@ void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
delete [] phase;
}
void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalty wetting strength W
{
// ONLY initialize grey nodes
// Key input parameters:
// 1. GreySolidLabels
// labels for grey nodes
// 2. GreySolidAffinity
// affinity ranges [-1,1]
// oil-wet > 0
// water-wet < 0
// neutral = 0
double *SolidPotential_host = new double [Nx*Ny*Nz];
double *GreySolidGrad_host = new double [3*Np];
// ranges [-1,1]
// water-wet > 0
// 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++){
@ -322,134 +336,48 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
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++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
Sn = SnList[idx];
Sw = SwList[idx];
idx = NLABELS;
//Mask->id[n] = 0; // set mask to zero since this is an immobile component
}
}
SolidPotential_host[n] = AFFINITY;
}
}
}
// Calculate grey-solid color-gradient
double *Dst;
Dst = new double [3*3*3];
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
}
}
}
double w_face = 1.f;
double w_edge = 0.5;
double w_corner = 0.f;
//local
Dst[13] = 0.f;
//faces
Dst[4] = w_face;
Dst[10] = w_face;
Dst[12] = w_face;
Dst[14] = w_face;
Dst[16] = w_face;
Dst[22] = w_face;
// corners
Dst[0] = w_corner;
Dst[2] = w_corner;
Dst[6] = w_corner;
Dst[8] = w_corner;
Dst[18] = w_corner;
Dst[20] = w_corner;
Dst[24] = w_corner;
Dst[26] = w_corner;
// edges
Dst[1] = w_edge;
Dst[3] = w_edge;
Dst[5] = w_edge;
Dst[7] = w_edge;
Dst[9] = w_edge;
Dst[11] = w_edge;
Dst[15] = w_edge;
Dst[17] = w_edge;
Dst[19] = w_edge;
Dst[21] = w_edge;
Dst[23] = w_edge;
Dst[25] = w_edge;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
int idx = Map(i,j,k);
if (!(idx < 0)){
double phi_x = 0.f;
double phi_y = 0.f;
double phi_z = 0.f;
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
double weight= Dst[index];
int idi=i+ii-1;
int idj=j+jj-1;
int idk=k+kk-1;
if (idi < 0) idi=0;
if (idj < 0) idj=0;
if (idk < 0) idk=0;
if (!(idi < Nx)) idi=Nx-1;
if (!(idj < Ny)) idj=Ny-1;
if (!(idk < Nz)) idk=Nz-1;
int nn = idk*Nx*Ny + idj*Nx + idi;
double vec_x = double(ii-1);
double vec_y = double(jj-1);
double vec_z = double(kk-1);
double GWNS=SolidPotential_host[nn];
phi_x += GWNS*weight*vec_x;
phi_y += GWNS*weight*vec_y;
phi_z += GWNS*weight*vec_z;
}
}
}
if (Averages->SDs(i,j,k)<2.0){
GreySolidGrad_host[idx+0*Np] = phi_x;
GreySolidGrad_host[idx+1*Np] = phi_y;
GreySolidGrad_host[idx+2*Np] = phi_z;
}
else{
GreySolidGrad_host[idx+0*Np] = 0.0;
GreySolidGrad_host[idx+1*Np] = 0.0;
GreySolidGrad_host[idx+2*Np] = 0.0;
}
}
GreySolidW_host[idx] = AFFINITY;
GreySn_host[idx] = Sn;
GreySw_host[idx] = Sw;
}
}
}
}
if (rank==0){
printf("Number of Grey-solid labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
AFFINITY=AffinityList[idx];
printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
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(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
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 [] SolidPotential_host;
delete [] GreySolidGrad_host;
delete [] Dst;
delete [] GreySolidW_host;
delete [] GreySn_host;
delete [] GreySw_host;
}
////----------------------------------------------------------------------------------------------------------//
@ -575,6 +503,71 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
delete [] Permeability;
}
//void ScaLBL_GreyscaleColorModel::AssignGreyscalePotential()
//{
// double *psi;//greyscale potential
// psi = new double[N];
//
// size_t NLABELS=0;
// signed char VALUE=0;
// double AFFINITY=0.f;
//
// auto LabelList = greyscaleColor_db->getVector<int>( "ComponentLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "ComponentAffinity" );
// NLABELS=LabelList.size();
//
// //first, copy over normal phase field
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// idx = NLABELS;
// }
// }
// // fluid labels are reserved
// if (VALUE == 1) AFFINITY=1.0;
// else if (VALUE == 2) AFFINITY=-1.0;
// psi[n] = AFFINITY;
// }
// }
// }
//
// //second, scale the phase field for grey nodes
// double Cap_Penalty=1.f;
// auto GreyLabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto PermeabilityList = greyscaleColor_db->getVector<double>( "PermeabilityList" );
// NLABELS=GreyLabelList.size();
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// Cap_Penalty=1.f;
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// if (VALUE == GreyLabelList[idx]){
// Cap_Penalty=alpha*W/sqrt(PermeabilityList[idx]/Dm->voxel_length/Dm->voxel_length);
// idx = NLABELS;
// }
// }
// //update greyscale potential
// psi[n] = psi[n]*Cap_Penalty;
// }
// }
// }
//
// ScaLBL_CopyToDevice(Psi, psi, N*sizeof(double));
// ScaLBL_Comm->Barrier();
// delete [] psi;
//}
void ScaLBL_GreyscaleColorModel::Create(){
/*
* This function creates the variables needed to run a LBM
@ -619,11 +612,15 @@ void ScaLBL_GreyscaleColorModel::Create(){
ScaLBL_AllocateDeviceMemory((void **) &Bq, 7*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nx*Ny*Nz);
//ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz);//greyscale potential
ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
//ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np);
//ScaLBL_AllocateDeviceMemory((void **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np);
//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);
//...........................................................................
@ -667,6 +664,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
AssignComponentLabels();//do open/black/grey nodes initialization
AssignGreySolidLabels();
AssignGreyPoroPermLabels();
//AssignGreyscalePotential();
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta,GreyPorosity);
ScaLBL_Comm->RegularLayout(Map,Porosity_dvc,Averages->Porosity);//porosity doesn't change over time
}
@ -931,9 +929,11 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Read for Aq, Bq happens in this routine (requires communication)
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_Update_GreyscalePotential(dvcMap,Phi,Psi,Porosity_dvc,Permeability_dvc,alpha,W,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
//ScaLBL_Update_GreyscalePotential(dvcMap,Phi,Psi,Porosity_dvc,Permeability_dvc,alpha,W,0,ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
@ -943,10 +943,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
// Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//Model-1&4
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
////Model-2&3
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
@ -968,10 +972,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
//Model-1&4
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
//Model-1&4
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
////Model-2&3
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
@ -983,9 +991,11 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Compute the Phase indicator field
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_Update_GreyscalePotential(dvcMap,Phi,Psi,Porosity_dvc,Permeability_dvc,alpha,W,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
//ScaLBL_Update_GreyscalePotential(dvcMap,Phi,Psi,Porosity_dvc,Permeability_dvc,alpha,W,0,ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
@ -995,10 +1005,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//Model-1&4
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
////Model-2&3
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
@ -1020,10 +1034,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
//Model-1&4
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
//Model-1&4
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
////Model-2&3
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
@ -1575,6 +1593,13 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE);
//ScaLBL_CopyToHost(PhaseField.data(), Psi, sizeof(double)*N);
//FILE *PSIFILE;
//sprintf(LocalRankFilename,"Psi.%05i.raw",rank);
//PSIFILE = fopen(LocalRankFilename,"wb");
//fwrite(PhaseField.data(),8,N,PSIFILE);
//fclose(PSIFILE);
ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
FILE *AFILE;
sprintf(LocalRankFilename,"A.%05i.raw",rank);
@ -1631,26 +1656,26 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
fwrite(PhaseField.data(),8,N,PERM_FILE);
fclose(PERM_FILE);
ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[0],PhaseField);
FILE *GreySG_X_FILE;
sprintf(LocalRankFilename,"GreySolidGrad_X.%05i.raw",rank);
GreySG_X_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,GreySG_X_FILE);
fclose(GreySG_X_FILE);
//ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[0],PhaseField);
//FILE *GreySG_X_FILE;
//sprintf(LocalRankFilename,"GreySolidGrad_X.%05i.raw",rank);
//GreySG_X_FILE = fopen(LocalRankFilename,"wb");
//fwrite(PhaseField.data(),8,N,GreySG_X_FILE);
//fclose(GreySG_X_FILE);
ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[Np],PhaseField);
FILE *GreySG_Y_FILE;
sprintf(LocalRankFilename,"GreySolidGrad_Y.%05i.raw",rank);
GreySG_Y_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,GreySG_Y_FILE);
fclose(GreySG_Y_FILE);
//ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[Np],PhaseField);
//FILE *GreySG_Y_FILE;
//sprintf(LocalRankFilename,"GreySolidGrad_Y.%05i.raw",rank);
//GreySG_Y_FILE = fopen(LocalRankFilename,"wb");
//fwrite(PhaseField.data(),8,N,GreySG_Y_FILE);
//fclose(GreySG_Y_FILE);
ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[2*Np],PhaseField);
FILE *GreySG_Z_FILE;
sprintf(LocalRankFilename,"GreySolidGrad_Z.%05i.raw",rank);
GreySG_Z_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,GreySG_Z_FILE);
fclose(GreySG_Z_FILE);
//ScaLBL_Comm->RegularLayout(Map,&GreySolidGrad[2*Np],PhaseField);
//FILE *GreySG_Z_FILE;
//sprintf(LocalRankFilename,"GreySolidGrad_Z.%05i.raw",rank);
//GreySG_Z_FILE = fopen(LocalRankFilename,"wb");
//fwrite(PhaseField.data(),8,N,GreySG_Z_FILE);
//fclose(GreySG_Z_FILE);
/* ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField);
FILE *CGX_FILE;
@ -1984,6 +2009,169 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
// delete [] Dst;
//}
//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
//{
// // ONLY initialize grey nodes
// // Key input parameters:
// // 1. GreySolidLabels
// // labels for grey nodes
// // 2. GreySolidAffinity
// // affinity ranges [-1,1]
// // oil-wet > 0
// // water-wet < 0
// // neutral = 0
// double *SolidPotential_host = new double [Nx*Ny*Nz];
// double *GreySolidGrad_host = new double [3*Np];
//
// size_t NLABELS=0;
// signed char VALUE=0;
// double AFFINITY=0.f;
//
// auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
//
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
// }
//
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// VALUE=id[n];
// AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// idx = NLABELS;
// //Mask->id[n] = 0; // set mask to zero since this is an immobile component
// }
// }
// SolidPotential_host[n] = AFFINITY;
// }
// }
// }
//
// // Calculate grey-solid color-gradient
// double *Dst;
// Dst = new double [3*3*3];
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
// int index = kk*9+jj*3+ii;
// Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
// }
// }
// }
// double w_face = 1.f;
// double w_edge = 0.5;
// double w_corner = 0.f;
// //local
// Dst[13] = 0.f;
// //faces
// Dst[4] = w_face;
// Dst[10] = w_face;
// Dst[12] = w_face;
// Dst[14] = w_face;
// Dst[16] = w_face;
// Dst[22] = w_face;
// // corners
// Dst[0] = w_corner;
// Dst[2] = w_corner;
// Dst[6] = w_corner;
// Dst[8] = w_corner;
// Dst[18] = w_corner;
// Dst[20] = w_corner;
// Dst[24] = w_corner;
// Dst[26] = w_corner;
// // edges
// Dst[1] = w_edge;
// Dst[3] = w_edge;
// Dst[5] = w_edge;
// Dst[7] = w_edge;
// Dst[9] = w_edge;
// Dst[11] = w_edge;
// Dst[15] = w_edge;
// Dst[17] = w_edge;
// Dst[19] = w_edge;
// Dst[21] = w_edge;
// Dst[23] = w_edge;
// Dst[25] = w_edge;
//
// for (int k=1; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
// int idx=Map(i,j,k);
// if (!(idx < 0)){
// double phi_x = 0.f;
// double phi_y = 0.f;
// double phi_z = 0.f;
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
//
// int index = kk*9+jj*3+ii;
// double weight= Dst[index];
//
// int idi=i+ii-1;
// int idj=j+jj-1;
// int idk=k+kk-1;
//
// if (idi < 0) idi=0;
// if (idj < 0) idj=0;
// if (idk < 0) idk=0;
// if (!(idi < Nx)) idi=Nx-1;
// if (!(idj < Ny)) idj=Ny-1;
// if (!(idk < Nz)) idk=Nz-1;
//
// int nn = idk*Nx*Ny + idj*Nx + idi;
// double vec_x = double(ii-1);
// double vec_y = double(jj-1);
// double vec_z = double(kk-1);
// double GWNS=SolidPotential_host[nn];
// phi_x += GWNS*weight*vec_x;
// phi_y += GWNS*weight*vec_y;
// phi_z += GWNS*weight*vec_z;
// }
// }
// }
// if (Averages->SDs(i,j,k)<2.0){
// GreySolidGrad_host[idx+0*Np] = phi_x;
// GreySolidGrad_host[idx+1*Np] = phi_y;
// GreySolidGrad_host[idx+2*Np] = phi_z;
// }
// else{
// GreySolidGrad_host[idx+0*Np] = 0.0;
// GreySolidGrad_host[idx+1*Np] = 0.0;
// GreySolidGrad_host[idx+2*Np] = 0.0;
// }
// }
// }
// }
// }
//
//
// if (rank==0){
// printf("Number of Grey-solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// VALUE=LabelList[idx];
// AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
//
// ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
// ScaLBL_Comm->Barrier();
// delete [] SolidPotential_host;
// delete [] GreySolidGrad_host;
// delete [] Dst;
//}
//--------- This is another old version of calculating greyscale-solid color-gradient modification-------//
// **not working effectively, to be deprecated
//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()

11
models/GreyscaleColorModel.h
View File

@ -39,6 +39,8 @@ public:
double Fx,Fy,Fz,flux;
double din,dout,inletA,inletB,outletA,outletB;
double GreyPorosity;
bool RecoloringOff;//recoloring can be turn off for grey nodes if this is true
//double W;//wetting strength paramter for capillary pressure penalty for grey nodes
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
@ -48,7 +50,7 @@ public:
std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular;
std::shared_ptr<GreyPhaseAnalysis> Averages;
std::shared_ptr<GreyPhaseAnalysis> Averages;
// input database
std::shared_ptr<Database> db;
@ -64,12 +66,16 @@ public:
double *fq, *Aq, *Bq;
double *Den, *Phi;
//double *GreySolidPhi; //Model 2 & 3
double *GreySolidGrad;//Model 1 & 4
//double *GreySolidGrad;//Model 1 & 4
double *GreySolidW;
double *GreySn;
double *GreySw;
//double *ColorGrad;
double *Velocity;
double *Pressure;
double *Porosity_dvc;
double *Permeability_dvc;
//double *Psi;
private:
Utilities::MPI comm;
@ -86,6 +92,7 @@ private:
void AssignComponentLabels();
void AssignGreySolidLabels();
void AssignGreyPoroPermLabels();
//void AssignGreyscalePotential();
void ImageInit(std::string filename);
double MorphInit(const double beta, const double morph_delta);
double SeedPhaseField(const double seed_water_in_oil);

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){

82
models/IonModel.cpp
View File

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

4
models/IonModel.h
View File

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

67
models/PoissonSolver.cpp
View File

@ -5,19 +5,35 @@
#include "analysis/distance.h"
#include "common/ReadMicroCT.h"
ScaLBL_Poisson::ScaLBL_Poisson(int RANK, int NP, const Utilities::MPI& COMM):
rank(RANK), nprocs(NP),timestep(0),timestepMax(0),tau(0),k2_inv(0),tolerance(0),h(0),
epsilon0(0),epsilon0_LB(0),epsilonR(0),epsilon_LB(0),Vin(0),Vout(0),Nx(0),Ny(0),Nz(0),N(0),Np(0),analysis_interval(0),
chargeDen_dummy(0),WriteLog(0),nprocx(0),nprocy(0),nprocz(0),
BoundaryConditionInlet(0),BoundaryConditionOutlet(0),BoundaryConditionSolid(0),Lx(0),Ly(0),Lz(0),
Vin0(0),freqIn(0),t0_In(0),Vin_Type(0),Vout0(0),freqOut(0),t0_Out(0),Vout_Type(0),
TestPeriodic(0),TestPeriodicTime(0),TestPeriodicTimeConv(0),TestPeriodicSaveInterval(0),
comm(COMM)
static inline bool fileExists( const std::string &filename )
{
std::ifstream ifile( filename.c_str() );
return ifile.good();
}
ScaLBL_Poisson::~ScaLBL_Poisson(){
ScaLBL_Poisson::ScaLBL_Poisson(int RANK, int NP, const Utilities::MPI& COMM):
rank(RANK), TIMELOG(nullptr), nprocs(NP),timestep(0),timestepMax(0),tau(0),k2_inv(0),tolerance(0),h(0),
epsilon0(0),epsilon0_LB(0),epsilonR(0),epsilon_LB(0),Vin(0),Vout(0),Nx(0),Ny(0),Nz(0),N(0),Np(0),analysis_interval(0),
chargeDen_dummy(0),WriteLog(0),nprocx(0),nprocy(0),nprocz(0),
BoundaryConditionInlet(0),BoundaryConditionOutlet(0),BoundaryConditionSolid(0),Lx(0),Ly(0),Lz(0),
Vin0(0),freqIn(0),t0_In(0),Vin_Type(0),Vout0(0),freqOut(0),t0_Out(0),Vout_Type(0),
TestPeriodic(0),TestPeriodicTime(0),TestPeriodicTimeConv(0),TestPeriodicSaveInterval(0),
comm(COMM)
{
if ( rank == 0 ) {
bool WriteHeader = !fileExists( "PoissonSolver_Convergence.csv" );
TIMELOG = fopen("PoissonSolver_Convergence.csv","a+");
if (WriteHeader)
fprintf(TIMELOG,"Timestep Error\n");
}
}
ScaLBL_Poisson::~ScaLBL_Poisson()
{
if ( TIMELOG )
fclose( TIMELOG );
}
void ScaLBL_Poisson::ReadParams(string filename){
@ -98,7 +114,6 @@ void ScaLBL_Poisson::ReadParams(string filename){
h = domain_db->getScalar<double>( "voxel_length" );
}
//Re-calcualte model parameters if user updates input
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity
@ -218,20 +233,19 @@ void ScaLBL_Poisson::ReadInput(){
void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid)
{
size_t NLABELS=0;
signed char VALUE=0;
double AFFINITY=0.f;
auto LabelList = electric_db->getVector<int>( "SolidLabels" );
auto AffinityList = electric_db->getVector<double>( "SolidValues" );
NLABELS=LabelList.size();
size_t NLABELS = LabelList.size();
if (NLABELS != AffinityList.size()){
ERROR("Error: LB-Poisson Solver: SolidLabels and SolidValues must be the same length! \n");
}
double label_count[NLABELS];
double label_count_global[NLABELS];
std::vector<double> label_count( NLABELS, 0.0 );
std::vector<double> label_count_global( NLABELS, 0.0 );
// Assign the labels
for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
@ -596,26 +610,11 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
}
void ScaLBL_Poisson::getConvergenceLog(int timestep,double error){
if (rank==0){
bool WriteHeader=false;
TIMELOG = fopen("PoissonSolver_Convergence.csv","r");
if (TIMELOG != NULL)
fclose(TIMELOG);
else
WriteHeader=true;
TIMELOG = fopen("PoissonSolver_Convergence.csv","a+");
if (WriteHeader)
{
fprintf(TIMELOG,"Timestep Error\n");
fprintf(TIMELOG,"%i %.5g\n",timestep,error);
fflush(TIMELOG);
}
else {
fprintf(TIMELOG,"%i %.5g\n",timestep,error);
fflush(TIMELOG);
}
void ScaLBL_Poisson::getConvergenceLog(int timestep,double error){
if ( rank == 0 ) {
fprintf(TIMELOG,"%i %.5g\n",timestep,error);
fflush(TIMELOG);
}
}

3
models/PoissonSolver.h
View File

@ -89,12 +89,13 @@ public:
private:
Utilities::MPI comm;
FILE *TIMELOG;
// filenames
char LocalRankString[8];
char LocalRankFilename[40];
char LocalRestartFile[40];
char OutputFilename[200];
FILE *TIMELOG;
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);

220
models/StokesModel.cpp
View File

@ -8,6 +8,7 @@
ScaLBL_StokesModel::ScaLBL_StokesModel(int RANK, int NP, const Utilities::MPI& COMM):
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),mu(0),h(0),nu_phys(0),rho_phys(0),rho0(0),den_scale(0),time_conv(0),tolerance(0),
epsilon0(0),epsilon0_LB(0),epsilonR(0),epsilon_LB(0),UseSlippingVelBC(0),
Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{
@ -38,6 +39,12 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
tolerance = 1.0e-8;
Fx = Fy = 0.0;
Fz = 1.0e-5;
//Stokes solver also needs the following parameters for slipping velocity BC
epsilon0 = 8.85e-12;//electric permittivity of vaccum; unit:[C/(V*m)]
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilonR = 78.4;//default dielectric constant of water
epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity
UseSlippingVelBC = false;
//--------------------------------------------------------------------------//
// Read domain parameters
@ -52,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" );
@ -85,22 +93,29 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
if (stokes_db->keyExists( "flux" )){
flux = stokes_db->getScalar<double>( "flux" );
}
if (stokes_db->keyExists( "UseElectroosmoticVelocityBC" )){
UseSlippingVelBC = stokes_db->getScalar<bool>( "UseElectroosmoticVelocityBC" );
}
if (stokes_db->keyExists( "epsilonR" )){
epsilonR = stokes_db->getScalar<double>( "epsilonR" );
}
// Re-calculate model parameters due to parameter read
mu=(tau-0.5)/3.0;
time_conv = (h*h*1.0e-12)*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity
}
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" );
stokes_db = db->getDatabase( "Stokes" );
//---------------------- Default model parameters --------------------------//
rho_phys = 1000.0; //by default use water density; unit [kg/m^3]
@ -114,6 +129,12 @@ void ScaLBL_StokesModel::ReadParams(string filename){
tolerance = 1.0e-8;
Fx = Fy = 0.0;
Fz = 1.0e-5;
//Stokes solver also needs the following parameters for slipping velocity BC
epsilon0 = 8.85e-12;//electric permittivity of vaccum; unit:[C/(V*m)]
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilonR = 78.4;//default dielectric constant of water
epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity
UseSlippingVelBC = false;
//--------------------------------------------------------------------------//
// Read domain parameters
@ -161,12 +182,19 @@ void ScaLBL_StokesModel::ReadParams(string filename){
if (stokes_db->keyExists( "flux" )){
flux = stokes_db->getScalar<double>( "flux" );
}
if (stokes_db->keyExists( "UseElectroosmoticVelocityBC" )){
UseSlippingVelBC = stokes_db->getScalar<bool>( "UseElectroosmoticVelocityBC" );
}
if (stokes_db->keyExists( "epsilonR" )){
epsilonR = stokes_db->getScalar<double>( "epsilonR" );
}
// Re-calculate model parameters due to parameter read
mu=(tau-0.5)/3.0;
time_conv = (h*h*1.0e-12)*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity
}
void ScaLBL_StokesModel::SetDomain(){
@ -258,6 +286,161 @@ void ScaLBL_StokesModel::ReadInput(){
if (rank == 0) cout << " Domain set." << endl;
}
void ScaLBL_StokesModel::AssignZetaPotentialSolid(double *zeta_potential_solid)
{
size_t NLABELS=0;
signed char VALUE=0;
double AFFINITY=0.f;
auto LabelList = stokes_db->getVector<int>( "SolidLabels" );
auto AffinityList = stokes_db->getVector<double>( "ZetaPotentialSolidList" );
NLABELS=LabelList.size();
if (NLABELS != AffinityList.size()){
ERROR("Error: LB Stokes Solver: SolidLabels and ZetaPotentialSolidList must be the same length! \n");
}
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;
// Assign the labels
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
VALUE=Mask->id[n];
AFFINITY=0.f;
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];//no need to convert unit for zeta potential (i.e. volt)
label_count[idx] += 1.0;
idx = NLABELS;
}
}
zeta_potential_solid[n] = AFFINITY;
}
}
}
for (size_t idx=0; idx<NLABELS; idx++)
label_count_global[idx]=Dm->Comm.sumReduce( label_count[idx]);
if (rank==0){
printf("LB Stokes Solver: number of solid labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
AFFINITY=AffinityList[idx];
double volume_fraction = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
printf(" label=%d, zeta potential=%.3g [V], volume fraction=%.2g\n",VALUE,AFFINITY,volume_fraction);
}
}
}
void ScaLBL_StokesModel::AssignSolidGrad(double *solid_grad)
{
double *Dst;
Dst = new double [3*3*3];
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
}
}
}
//implement a D3Q19 lattice
double w_face = 1.0/18.0;
double w_edge = 0.5*w_face;
double w_corner = 0.0;
//local
Dst[13] = 0.f;
//faces
Dst[4] = w_face;
Dst[10] = w_face;
Dst[12] = w_face;
Dst[14] = w_face;
Dst[16] = w_face;
Dst[22] = w_face;
// corners
Dst[0] = w_corner;
Dst[2] = w_corner;
Dst[6] = w_corner;
Dst[8] = w_corner;
Dst[18] = w_corner;
Dst[20] = w_corner;
Dst[24] = w_corner;
Dst[26] = w_corner;
// edges
Dst[1] = w_edge;
Dst[3] = w_edge;
Dst[5] = w_edge;
Dst[7] = w_edge;
Dst[9] = w_edge;
Dst[11] = w_edge;
Dst[15] = w_edge;
Dst[17] = w_edge;
Dst[19] = w_edge;
Dst[21] = w_edge;
Dst[23] = w_edge;
Dst[25] = w_edge;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
if (!(idx < 0)){
double phi_x = 0.f;
double phi_y = 0.f;
double phi_z = 0.f;
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
double weight= Dst[index];
int idi=i+ii-1;
int idj=j+jj-1;
int idk=k+kk-1;
if (idi < 0) idi=0;
if (idj < 0) idj=0;
if (idk < 0) idk=0;
if (!(idi < Nx)) idi=Nx-1;
if (!(idj < Ny)) idj=Ny-1;
if (!(idk < Nz)) idk=Nz-1;
int nn = idk*Nx*Ny + idj*Nx + idi;
double vec_x = double(ii-1);
double vec_y = double(jj-1);
double vec_z = double(kk-1);
double GWNS = double(Mask->id[nn]);
//Since the solid unit normal vector is wanted, treat
//wet node as 0.0 and solid node as 1.0
GWNS = (GWNS>0.0) ? 0.0:1.0;
phi_x += GWNS*weight*vec_x;
phi_y += GWNS*weight*vec_y;
phi_z += GWNS*weight*vec_z;
}
}
}
//solid_grad normalization
double phi_mag=sqrt(phi_x*phi_x+phi_y*phi_y+phi_z*phi_z);
if (phi_mag==0.0) phi_mag=1.0;
solid_grad[idx+0*Np] = phi_x/phi_mag;
solid_grad[idx+1*Np] = phi_y/phi_mag;
solid_grad[idx+2*Np] = phi_z/phi_mag;
}
}
}
}
}
void ScaLBL_StokesModel::Create(){
/*
* This function creates the variables needed to run a LBM
@ -301,6 +484,26 @@ void ScaLBL_StokesModel::Create(){
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
comm.barrier();
if (UseSlippingVelBC==true){
ScaLBL_Comm->SetupBounceBackList(Map, Mask->id.data(), Np,1);
comm.barrier();
//For slipping velocity BC, need zeta potential and solid unit normal vector
ScaLBL_AllocateDeviceMemory((void **) &ZetaPotentialSolid, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &SolidGrad, sizeof(double)*3*Np); //unit normal vector of solid nodes
double *ZetaPotentialSolid_host;
ZetaPotentialSolid_host = new double[Nx*Ny*Nz];
AssignZetaPotentialSolid(ZetaPotentialSolid_host);
double *SolidGrad_host;
SolidGrad_host = new double[3*Np];
AssignSolidGrad(SolidGrad_host);
ScaLBL_CopyToDevice(ZetaPotentialSolid, ZetaPotentialSolid_host, Nx*Ny*Nz*sizeof(double));
ScaLBL_CopyToDevice(SolidGrad, SolidGrad_host, 3*Np*sizeof(double));
ScaLBL_Comm->Barrier();
delete [] ZetaPotentialSolid_host;
delete [] SolidGrad_host;
}
}
void ScaLBL_StokesModel::Initialize(){
@ -324,6 +527,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
timestep = 0;
while (timestep < timestepMax) {
//************************************************************************/
//**************ODD TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
@ -344,8 +548,14 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
}
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
0, ScaLBL_Comm->LastExterior(), Np);
if (UseSlippingVelBC==true){
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(fq, ZetaPotentialSolid, ElectricField, SolidGrad,
epsilon_LB, 1.0/rlx_setA, rho0, den_scale, h, time_conv);
}
ScaLBL_Comm->Barrier(); comm.barrier();
//**************EVEN TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
@ -366,6 +576,10 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
}
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
0, ScaLBL_Comm->LastExterior(), Np);
if (UseSlippingVelBC==true){
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(fq, ZetaPotentialSolid, ElectricField, SolidGrad,
epsilon_LB, 1.0/rlx_setA, rho0, den_scale, h, time_conv);
}
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
}

6
models/StokesModel.h
View File

@ -51,6 +51,8 @@ public:
double time_conv;
double h;//image resolution
double den_scale;//scale factor for density
double epsilon0,epsilon0_LB,epsilonR,epsilon_LB;//Stokes solver also needs this for slipping velocity BC
bool UseSlippingVelBC;
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
@ -70,6 +72,8 @@ public:
double *fq;
double *Velocity;
double *Pressure;
double *ZetaPotentialSolid;
double *SolidGrad;
//Minkowski Morphology;
DoubleArray Velocity_x;
@ -88,5 +92,7 @@ private:
void LoadParams(std::shared_ptr<Database> db0);
void Velocity_LB_to_Phys(DoubleArray &Vel_reg);
vector<double> computeElectricForceAvg(double *ChargeDensity, double *ElectricField);
void AssignSolidGrad(double *solid_grad);
void AssignZetaPotentialSolid(double *zeta_potential_solid);
};
#endif

34
sample_scripts/configure_lonestar Executable file
View File

@ -0,0 +1,34 @@
# load the module for cmake
#module load cmake
#source /gpfs/gpfs_stage1/b6p315aa/setup/setup-mpi.sh
module load cmake gcc
module load cuda
export HDF5_DIR=$HOME/local/hdf5/1.8.12/
export SILO_DIR=$HOME/local/silo/4.10.2/
export NETCDF_DIR=$HOME/local/netcdf/4.6.1
# configure
rm -rf CMake*
cmake \
-D CMAKE_BUILD_TYPE:STRING=Release \
-D CMAKE_C_COMPILER:PATH=mpicc \
-D CMAKE_CXX_COMPILER:PATH=mpicxx \
-D CMAKE_CXX_STANDARD=14 \
-D USE_CUDA=1 \
-D CMAKE_CUDA_FLAGS="-arch sm_70 -Xptxas=-v -Xptxas -dlcm=cg -lineinfo" \
-D CMAKE_CUDA_HOST_COMPILER="/opt/apps/gcc/7.3.0/bin/gcc" \
-D USE_HDF5=1 \
-D HDF5_DIRECTORY="$HDF5_DIR" \
-D HDF5_LIB="$HDF5_DIR/lib/libhdf5.a" \
-D USE_SILO=1 \
-D SILO_LIB="$SILO_DIR/lib/libsiloh5.a" \
-D SILO_DIRECTORY="$SILO_DIR" \
-D USE_NETCDF=0 \
-D NETCDF_DIRECTORY="$NETCDF_DIR" \
-D USE_DOXYGEN:BOOL=false \
-D USE_TIMER=0 \
~/src/LBPM
make VERBOSE=1 -j1 && make install

2
sample_scripts/configure_summit
View File

@ -4,7 +4,7 @@
#source /gpfs/gpfs_stage1/b6p315aa/setup/setup-mpi.sh
module load cmake gcc
module load cuda
#/ccs/proj/csc380/mcclurej
export HDF5_DIR=/ccs/proj/csc380/mcclurej/install/hdf5/1.8.12/
export SILO_DIR=/ccs/proj/csc380/mcclurej/install/silo/4.10.2/
export NETCDF_DIR=/ccs/proj/geo136/install/netcdf/4.6.1

2
tests/CMakeLists.txt
View File

@ -62,7 +62,7 @@ ADD_LBPM_TEST( TestMap )
ADD_LBPM_TEST( TestWideHalo )
ADD_LBPM_TEST( TestColorGradDFH )
ADD_LBPM_TEST( TestBubbleDFH ../example/Bubble/input.db)
ADD_LBPM_TEST( testGlobalMassFreeLee ../example/Bubble/input.db)
#ADD_LBPM_TEST( testGlobalMassFreeLee ../example/Bubble/input.db)
#ADD_LBPM_TEST( TestColorMassBounceback ../example/Bubble/input.db)
ADD_LBPM_TEST( TestPressVel ../example/Bubble/input.db)
ADD_LBPM_TEST( TestPoiseuille ../example/Piston/poiseuille.db)

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;

226
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,82 +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);
}
} //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();
ColorModel.WriteDebug();
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 );
// ****************************************************

3
tests/lbpm_freelee_simulator.cpp
View File

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