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

merging conflicts

Browse Source
This commit is contained in:
James McClure
2021-12-09 13:51:35 -05:00
parent 8800066c06 ee93851281
commit 51ada19b06
10 changed files with 3288 additions and 4180 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4839
common/ScaLBL.cpp
View File

File diff suppressed because it is too large Load Diff

891
common/ScaLBL.h
View File

File diff suppressed because it is too large Load Diff

32
cpu/D3Q7BC.cpp
View File

@@ -35,13 +35,31 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValue,
}
}
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_Solid_DirichletAndNeumann_D3Q7(double *dist,double *BoundaryValue,int* BoundaryLabel,int *BounceBackDist_list,int *BounceBackSolid_list,int N){
int idx;
int iq,ib;
double value_b,value_b_label,value_q;
for (idx=0; idx<N; idx++){
iq = BounceBackDist_list[idx];
ib = BounceBackSolid_list[idx];
value_b = BoundaryValue[ib];//get boundary value from a solid site
value_b_label = BoundaryLabel[ib];//get boundary label (i.e. type of BC) from a solid site
value_q = dist[iq];
if (value_b_label==1){//Dirichlet BC
dist[iq] = -1.0*value_q + value_b*0.25;//NOTE 0.25 is the speed of sound for D3Q7 lattice
}
if (value_b_label==2){//Neumann BC
dist[iq] = value_q + value_b;
}
}
}
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;

31
cuda/D3Q7BC.cu
View File

@@ -38,6 +38,28 @@ __global__ void dvc_ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValu
}
}
__global__ void dvc_ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist, double *BoundaryValue,int *BoundaryLabel, int *BounceBackDist_list, int *BounceBackSolid_list, int count)
{
int idx;
int iq,ib;
double value_b,value_b_label,value_q;
idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
iq = BounceBackDist_list[idx];
ib = BounceBackSolid_list[idx];
value_b = BoundaryValue[ib];//get boundary value from a solid site
value_b_label = BoundaryLabel[ib];//get boundary label (i.e. type of BC) from a solid site
value_q = dist[iq];
if (value_b_label==1){//Dirichlet BC
dist[iq] = -1.0*value_q + value_b*0.25;//NOTE 0.25 is the speed of sound for D3Q7 lattice
}
if (value_b_label==2){//Neumann BC
dist[iq] = value_q + value_b;
}
}
}
__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,
@@ -733,6 +755,15 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValue, i
}
}
extern "C" void ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist, double *BoundaryValue,int *BoundaryLabel, int *BounceBackDist_list, int *BounceBackSolid_list, int count){
int GRID = count / 512 + 1;
dvc_ScaLBL_Solid_DirichletAndNeumann_D3Q7<<<GRID,512>>>(dist, BoundaryValue, BoundaryLabel, BounceBackDist_list, BounceBackSolid_list, count);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_Solid_DirichletAndNeumann_D3Q7 (kernel): %s \n",cudaGetErrorString(err));
}
}
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,

123
docs/source/examples/color/steadyState.rst
View File

@@ -1,5 +1,120 @@
*******************
Steady-state flow
*******************
********************************
Steady-state flow (color model)
********************************
In this example we simulate a steady-state flow with a constant driving force. This will enforce a periodic boundary condition
in all directions. While the driving force may be set in any direction, we will set it in the z-direction to be consistent
with the convention for pressure and velocity boundary conditions.
For the case considered in ``example/DiscPack`` we specify the following information in the input file
.. code:: c
Domain {
Filename = "discs_3x128x128.raw.morphdrain.raw"
ReadType = "8bit" // data type
N = 3, 128, 128 // size of original image
nproc = 1, 2, 2 // process grid
n = 3, 64, 64 // sub-domain size
voxel_length = 1.0 // voxel length (in microns)
ReadValues = 0, 1, 2 // labels within the original image
WriteValues = 0, 1, 2 // associated labels to be used by LBPM
BC = 0 // fully periodic BC
Sw = 0.35 // target saturation for morphological tools
}
Color {
capillary_number = -1e-5 // capillary number for the displacement, positive="oil injection"
timestepMax = 20000 // maximum timtestep
alpha = 0.01 // controls interfacial tension
rhoA = 1.0 // controls the density of fluid A
rhoB = 1.0 // controls the density of fluid B
tauA = 0.7 // controls the viscosity of fluid A
tauB = 0.7 // controls the viscosity of fluid B
F = 0, 0, 1e-5 // body force
WettingConvention = "SCAL"
ComponentLabels = 0 // image labels for solid voxels
ComponentAffinity = 0.9 // controls the wetting affinity for each label
Restart = false
}
Analysis {
analysis_interval = 1000 // logging interval for timelog.csv
subphase_analysis_interval = 500000 // loggging interval for subphase.csv
N_threads = 4 // number of analysis threads (GPU version only)
visualization_interval = 10000 // interval to write visualization files
restart_interval = 10000000 // interval to write restart file
restart_file = "Restart" // base name of restart file
}
Visualization {
format = "hdf5"
write_silo = true // write SILO databases with assigned variables
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
save_phase_field = true // save phase field within SILO database
save_pressure = true // save pressure field within SILO database
save_velocity = false // save velocity field within SILO database
}
FlowAdaptor {
}
Once this has been set, we launch ``lbpm_color_simulator`` in the same way as other parallel tools
.. code:: bash
mpirun -np 4 $LBPM_BIN/lbpm_color_simulator input.db
Successful output looks like the following
.. code:: bash
********************************************************
Running Color LBM
********************************************************
voxel length = 1.000000 micron
voxel length = 1.000000 micron
Input media: discs_3x128x128.raw.morphdrain.raw
Relabeling 3 values
oldvalue=0, newvalue =0
oldvalue=1, newvalue =1
oldvalue=2, newvalue =2
Dimensions of segmented image: 3 x 128 x 128
Reading 8-bit input data
Read segmented data from discs_3x128x128.raw.morphdrain.raw
Label=0, Count=11862
Label=1, Count=26430
Label=2, Count=10860
Distributing subdomains across 4 processors
Process grid: 1 x 2 x 2
Subdomain size: 3 x 64 x 64
Size of transition region: 0
Media porosity = 0.758667
Initialized solid phase -- Converting to Signed Distance function
Domain set.
Create ScaLBL_Communicator
Set up memory efficient layout, 9090 | 9120 | 21780
Allocating distributions
Setting up device map and neighbor list
Component labels: 1
label=0, affinity=-0.900000, volume fraction==0.417582
Initializing distributions
Initializing phase field
Affinities - rank 0:
Main: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Thread 1: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Thread 2: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Thread 3: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Thread 4: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Affinities - rank 0:
Main: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Thread 1: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Thread 2: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Thread 3: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
Thread 4: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
********************************************************
CPU time = 0.001501
Lattice update rate (per core)= 6.074861 MLUPS
Lattice update rate (per MPI process)= 6.074861 MLUPS
(flatten density field)
In this example we simulate a steady-state flow with a constant driving force.

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

@@ -225,4 +225,5 @@ Example Input File
InletLayers = 0, 0, 10 // specify 10 layers along the z-inlet
BC = 0 // boundary condition type (0 for periodic)
}
Visualization {
}

31
hip/D3Q7BC.cu
View File

@@ -37,6 +37,28 @@ __global__ void dvc_ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValu
}
}
__global__ void dvc_ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist, double *BoundaryValue,int *BoundaryLabel, int *BounceBackDist_list, int *BounceBackSolid_list, int count)
{
int idx;
int iq,ib;
double value_b,value_b_label,value_q;
idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
iq = BounceBackDist_list[idx];
ib = BounceBackSolid_list[idx];
value_b = BoundaryValue[ib];//get boundary value from a solid site
value_b_label = BoundaryLabel[ib];//get boundary label (i.e. type of BC) from a solid site
value_q = dist[iq];
if (value_b_label==1){//Dirichlet BC
dist[iq] = -1.0*value_q + value_b*0.25;//NOTE 0.25 is the speed of sound for D3Q7 lattice
}
if (value_b_label==2){//Neumann BC
dist[iq] = value_q + value_b;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np)
{
int idx,n;
@@ -409,6 +431,15 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValue, i
}
}
extern "C" void ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist, double *BoundaryValue,int *BoundaryLabel, int *BounceBackDist_list, int *BounceBackSolid_list, int count){
int GRID = count / 512 + 1;
dvc_ScaLBL_Solid_DirichletAndNeumann_D3Q7<<<GRID,512>>>(dist, BoundaryValue, BoundaryLabel, BounceBackDist_list, BounceBackSolid_list, count);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("hip error in ScaLBL_Solid_DirichletAndNeumann_D3Q7 (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);

38
models/ColorModel.cpp
View File

@@ -618,7 +618,6 @@ double ScaLBL_ColorModel::Run(int returntime) {
bool SET_CAPILLARY_NUMBER = false;
bool TRIGGER_FORCE_RESCALE = false;
double tolerance = 0.01;
auto WettingConvention = color_db->getWithDefault<std::string>( "WettingConvention", "none" );
auto current_db = db->cloneDatabase();
auto flow_db = db->getDatabase("FlowAdaptor");
int MIN_STEADY_TIMESTEPS =
@@ -646,7 +645,7 @@ double ScaLBL_ColorModel::Run(int returntime) {
if (analysis_db->keyExists("tolerance")) {
tolerance = analysis_db->getScalar<double>("tolerance");
}
runAnalysis analysis(current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular,
Map);
auto t1 = std::chrono::system_clock::now();
@@ -965,41 +964,6 @@ double ScaLBL_ColorModel::Run(int returntime) {
fprintf(kr_log_file, "%.5g\n", eff_pres);
fclose(kr_log_file);
if (WettingConvention == "SCAL"){
WriteHeader = false;
FILE *scal_log_file = fopen("SCAL.csv", "r");
if (scal_log_file != NULL)
fclose(scal_log_file);
else
WriteHeader = true;
scal_log_file = fopen("SCAL.csv", "a");
if (WriteHeader) {
fprintf(scal_log_file, "timesteps sat.water ");
fprintf(scal_log_file, "eff.perm.oil.upper.bound "
"eff.perm.water.upper.bound ");
fprintf(scal_log_file,
"eff.perm.oil.lower.bound "
"eff.perm.water.lower.bound ");
fprintf(scal_log_file, "eff.perm.oil.disconnected "
"eff.perm.water.disconnected ");
fprintf(scal_log_file,
"cap.pressure cap.pressure.connected "
"Ca eff.pressure\n");
}
fprintf(scal_log_file, "%i %.5g ", CURRENT_TIMESTEP,
current_saturation);
fprintf(scal_log_file, "%.5g %.5g ", kAeff_low, kBeff_low);
fprintf(scal_log_file, "%.5g %.5g ", kAeff_connected_low,
kBeff_connected_low);
fprintf(scal_log_file, "%.5g %.5g ", kAeff_disconnected,
kBeff_disconnected);
fprintf(scal_log_file, "%.5g %.5g %.5g ", pAB,
pAB_connected, Ca);
fprintf(scal_log_file, "%.5g\n", eff_pres);
fclose(scal_log_file);
}
printf(" Measured capillary number %f \n ", Ca);
}
if (SET_CAPILLARY_NUMBER) {

1447
models/PoissonSolver.cpp
View File

File diff suppressed because it is too large Load Diff

33
models/PoissonSolver.h
View File

@@ -45,22 +45,22 @@ public:
//bool Restart,pBC;
int timestep, timestepMax;
int analysis_interval;
int BoundaryConditionInlet;
int BoundaryConditionOutlet;
int BoundaryConditionSolid;
double tau;
double tolerance;
int BoundaryConditionInlet;
int BoundaryConditionOutlet;
vector<int> BoundaryConditionSolidList;
double tau;
double tolerance;
std::string tolerance_method;
double k2_inv;
double epsilon0, epsilon0_LB, epsilonR, epsilon_LB;
double Vin, Vout;
double chargeDen_dummy; //for debugging
bool WriteLog;
double Vin0, freqIn, t0_In, Vin_Type;
double Vout0, freqOut, t0_Out, Vout_Type;
bool TestPeriodic;
double TestPeriodicTime; //unit: [sec]
double TestPeriodicTimeConv; //unit [sec/lt]
double Vin0,freqIn,t0_In,PhaseShift_In;
double Vout0,freqOut,t0_Out,PhaseShift_Out;
bool TestPeriodic;
double TestPeriodicTime;//unit: [sec]
double TestPeriodicTimeConv; //unit [sec/lt]
double TestPeriodicSaveInterval; //unit [sec]
int Nx, Ny, Nz, N, Np;
@@ -86,7 +86,8 @@ public:
int *dvcMap;
//signed char *dvcID;
double *fq;
double *Psi;
double *Psi;
int *Psi_BCLabel;
double *ElectricField;
double *ChargeDensityDummy; // for debugging
double *ResidualError;
@@ -103,8 +104,8 @@ private:
char OutputFilename[200];
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void AssignSolidBoundary(double *poisson_solid);
void LoadParams(std::shared_ptr<Database> db0);
void AssignSolidBoundary(double *poisson_solid, int *poisson_solid_label);
void Potential_Init(double *psi_init);
void ElectricField_LB_to_Phys(DoubleArray &Efield_reg);
void SolveElectricPotentialAAodd(int timestep_from_Study);
@@ -112,8 +113,8 @@ private:
//void SolveElectricField();
void SolvePoissonAAodd(double *ChargeDensity);
void SolvePoissonAAeven(double *ChargeDensity);
void getConvergenceLog(int timestep, double error);
double getBoundaryVoltagefromPeriodicBC(double V0, double freq, double t0,
int V_type, int time_step);
void getConvergenceLog(int timestep,double error);
double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int time_step);
};
#endif