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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James E McClure 2022-03-05 07:47:22 -05:00
commit 77eba8c48a
21 changed files with 162 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
analysis/morphology.cpp
View File

@ -611,7 +611,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
//***************************************************************************************
double MorphDrain(DoubleArray &SignDist, signed char *id,
std::shared_ptr<Domain> Dm, double VoidFraction) {
std::shared_ptr<Domain> Dm, double VoidFraction, double InitialRadius) {
// SignDist is the distance to the object that you want to constaing the morphological opening
// VoidFraction is the the empty space where the object inst
// id is a labeled map
@ -688,6 +688,11 @@ double MorphDrain(DoubleArray &SignDist, signed char *id,
double deltaR = 0.05; // amount to change the radius in voxel units
double Rcrit_old = maxdistGlobal;
double Rcrit_new = maxdistGlobal;
if (InitialRadius < maxdistGlobal){
Rcrit_old = InitialRadius;
Rcrit_new = InitialRadius;
}
//if (argc>2){
// Rcrit_new = strtod(argv[2],NULL);
// if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new);

2
analysis/morphology.h
View File

@ -7,7 +7,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
std::shared_ptr<Domain> Dm, double VoidFraction,
signed char ErodeLabel, signed char ReplaceLabel);
double MorphDrain(DoubleArray &SignDist, signed char *id,
std::shared_ptr<Domain> Dm, double VoidFraction);
std::shared_ptr<Domain> Dm, double VoidFraction, double InitialRadius);
double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
std::shared_ptr<Domain> Dm, double TargetVol,
double WallFactor);

13
common/ScaLBL.h
View File

@ -296,7 +296,7 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB, bool UseSlippingVelBC,
int start, int finish, int Np);
/**
@ -307,19 +307,22 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *di
* @param Psi -
* @param ElectricField - electric field
* @param tau - relaxation time
* @param epsilon_LB -
* @param epsilon_LB - dielectric constant of medium
* @param UseSlippingVelBC - flag indicating the use of Helmholtz-Smoluchowski slipping velocity equation when EDL is too small to resolve
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
*/
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau,
double epsilon_LB, int start, int finish, int Np);
double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np);
/**
* \brief Poisson-Boltzmann solver / solve electric potential based on AA odd access pattern for D3Q7
* @param neighborList - neighbors based on D3Q19 lattice structure
* @param Map - mapping between sparse and dense representations
* @param dist - D3Q7 distributions
* @param Psi -
* @param epsilon_LB - dielectric constant of medium
* @param UseSlippingVelBC - flag indicating the use of Helmholtz-Smoluchowski slipping velocity equation when EDL is too small to resolve
* @param start - lattice node to start loop
* @param finish - lattice node to finish loop
* @param Np - size of local sub-domain (derived from Domain structure)
@ -350,10 +353,10 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, in
// LBM Stokes Model (adapted from MRT model)
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np);
double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np);
double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np);
extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np);

16
cpu/Poisson.cpp
View File

@ -95,7 +95,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
double *dist, double *Den_charge,
double *Psi, double *ElectricField,
double tau, double epsilon_LB,
double tau, double epsilon_LB, bool UseSlippingVelBC,
int start, int finish, int Np) {
int n;
double psi; //electric potential
@ -109,8 +109,9 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
for (n = start; n < finish; n++) {
//Load data
rho_e = Den_charge[n];
rho_e = rho_e / epsilon_LB;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
idx = Map[n];
psi = Psi[idx];
@ -175,8 +176,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist,
double *Den_charge, double *Psi,
double *ElectricField, double tau,
double epsilon_LB, int start,
int finish, int Np) {
double epsilon_LB, bool UseSlippingVelBC,
int start, int finish, int Np) {
int n;
double psi; //electric potential
double Ex, Ey, Ez; //electric field
@ -188,8 +189,9 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist,
for (n = start; n < finish; n++) {
//Load data
rho_e = Den_charge[n];
rho_e = rho_e / epsilon_LB;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
idx = Map[n];
psi = Psi[idx];

31
cpu/Stokes.cpp
View File

@ -4,7 +4,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
double *dist, double *Velocity, double *ChargeDensity,
double *ElectricField, double rlx_setA, double rlx_setB, double Gx,
double Gy, double Gz, double rho0, double den_scale, double h,
double time_conv, int start, int finish, int Np) {
double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) {
double fq;
// conserved momemnts
double rho, jx, jy, jz;
@ -38,13 +38,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
//compute total body force, including input body force (Gx,Gy,Gz)
Fx =
Gx +
rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
// q=0
@ -479,7 +477,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(
int *neighborList, double *dist, double *Velocity, double *ChargeDensity,
double *ElectricField, double rlx_setA, double rlx_setB, double Gx,
double Gy, double Gz, double rho0, double den_scale, double h,
double time_conv, int start, int finish, int Np) {
double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) {
double fq;
// conserved momemnts
double rho, jx, jy, jz;
@ -513,12 +511,21 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(
Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
//compute total body force, including input body force (Gx,Gy,Gz)
Fx = Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
//Fx = Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
// den_scale; //the extra factors at the end necessarily convert unit from phys to LB
//Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
// den_scale;
//Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
// den_scale;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and body force induced by external efectric field is reduced to slipping velocity BC.
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
// q=0

22
cuda/Poisson.cu
View File

@ -104,7 +104,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, doub
}
}
__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
int n;
double psi;//electric potential
@ -122,8 +122,9 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
if (n<finish) {
//Load data
rho_e = Den_charge[n];
rho_e = rho_e/epsilon_LB;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
idx=Map[n];
psi = Psi[idx];
@ -184,7 +185,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
int n;
double psi;//electric potential
@ -201,8 +202,9 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *
if (n<finish) {
//Load data
rho_e = Den_charge[n];
rho_e = rho_e/epsilon_LB;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
idx=Map[n];
psi = Psi[idx];
@ -293,10 +295,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
@ -305,10 +307,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){

30
cuda/Stokes.cu
View File

@ -5,7 +5,7 @@
#define NBLOCKS 1024
#define NTHREADS 256
__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,bool UseSlippingVelBC,int start, int finish, int Np){
int n;
double fq;
@ -46,9 +46,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
Ey = ElectricField[n+1*Np];
Ez = ElectricField[n+2*Np];
//compute total body force, including input body force (Gx,Gy,Gz)
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
// q=0
fq = dist[n];
@ -510,7 +513,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
}
}
__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
int n;
double fq;
@ -550,9 +553,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
Ey = ElectricField[n+1*Np];
Ez = ElectricField[n+2*Np];
//compute total body force, including input body force (Gx,Gy,Gz)
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
// q=0
fq = dist[n];
@ -969,10 +975,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
}
}
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC,start,finish,Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
@ -981,10 +987,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC,start,finish,Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){

22
hip/Poisson.cu
View File

@ -104,7 +104,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, doub
}
}
__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
int n;
double psi;//electric potential
@ -122,8 +122,9 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
if (n<finish) {
//Load data
rho_e = Den_charge[n];
rho_e = rho_e/epsilon_LB;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
idx=Map[n];
psi = Psi[idx];
@ -184,7 +185,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
}
}
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
int n;
double psi;//electric potential
@ -201,8 +202,9 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *
if (n<finish) {
//Load data
rho_e = Den_charge[n];
rho_e = rho_e/epsilon_LB;
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
idx=Map[n];
psi = Psi[idx];
@ -293,10 +295,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
@ -305,10 +307,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){

30
hip/Stokes.cu
View File

@ -6,7 +6,7 @@
#define NTHREADS 256
__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,bool UseSlippingVelBC,int start, int finish, int Np){
int n;
double fq;
@ -47,9 +47,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
Ey = ElectricField[n+1*Np];
Ez = ElectricField[n+2*Np];
//compute total body force, including input body force (Gx,Gy,Gz)
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
// q=0
fq = dist[n];
@ -511,7 +514,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
}
}
__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC,int start, int finish, int Np){
int n;
double fq;
@ -551,9 +554,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
Ey = ElectricField[n+1*Np];
Ez = ElectricField[n+2*Np];
//compute total body force, including input body force (Gx,Gy,Gz)
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
den_scale;
// q=0
fq = dist[n];
@ -970,10 +976,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
}
}
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
@ -982,10 +988,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){

16
models/ColorModel.cpp
View File

@ -876,6 +876,14 @@ double ScaLBL_ColorModel::Run(int returntime) {
double vBd_x = Averages->gwd.Px / Mass_w;
double vBd_y = Averages->gwd.Py / Mass_w;
double vBd_z = Averages->gwd.Pz / Mass_w;
/* contribution from water films */
double water_film_flow_rate =
Averages->gwb.V * (Averages->giwn.Pwx * dir_x + Averages->giwn.Pwy * dir_y + Averages->giwn.Pwz * dir_z) /
Averages->gwb.M / Dm->Volume;
double not_water_film_flow_rate =
Averages->gnb.V * (Averages->giwn.Pnx * dir_x + Averages->giwn.Pny * dir_y + Averages->giwn.Pnz * dir_z) /
Averages->gnb.M / Dm->Volume;
double flow_rate_A_connected =
Vol_nc *
@ -912,6 +920,11 @@ double ScaLBL_ColorModel::Run(int returntime) {
double kAeff = h * h * muA * (flow_rate_A) / (force_mag);
double kBeff = h * h * muB * (flow_rate_B) / (force_mag);
/* not counting films */
double krnf = kAeff - h * h * muA * not_water_film_flow_rate / force_mag;
double krwf = kBeff - h * h * muB * water_film_flow_rate / force_mag;
// Saturation normalized effective permeability to account for decoupled phases and
// effective porosity.
double kAeff_low = (1.0 - current_saturation) * h * h *
@ -936,6 +949,8 @@ double ScaLBL_ColorModel::Run(int returntime) {
fprintf(kr_log_file, "timesteps sat.water ");
fprintf(kr_log_file, "eff.perm.oil.upper.bound "
"eff.perm.water.upper.bound ");
fprintf(kr_log_file, "eff.perm.oil.film "
"eff.perm.water.film ");
fprintf(kr_log_file, "eff.perm.oil.lower.bound "
"eff.perm.water.lower.bound ");
fprintf(kr_log_file,
@ -953,6 +968,7 @@ double ScaLBL_ColorModel::Run(int returntime) {
fprintf(kr_log_file, "%i %.5g ", CURRENT_TIMESTEP,
current_saturation);
fprintf(kr_log_file, "%.5g %.5g ", kAeff, kBeff);
fprintf(kr_log_file, "%.5g %.5g ", krnf, krwf);
fprintf(kr_log_file, "%.5g %.5g ", kAeff_low, kBeff_low);
fprintf(kr_log_file, "%.5g %.5g ", kAeff_connected,
kBeff_connected);

3
models/MultiPhysController.cpp
View File

@ -4,7 +4,7 @@ ScaLBL_Multiphys_Controller::ScaLBL_Multiphys_Controller(
int RANK, int NP, const Utilities::MPI &COMM)
: rank(RANK), nprocs(NP), Restart(0), timestepMax(0), num_iter_Stokes(0),
num_iter_Ion(0), analysis_interval(0), visualization_interval(0),
tolerance(0), time_conv_max(0), comm(COMM) {}
tolerance(0), time_conv_max(0), time_conv_MainLoop(0), comm(COMM) {}
ScaLBL_Multiphys_Controller::~ScaLBL_Multiphys_Controller() {}
void ScaLBL_Multiphys_Controller::ReadParams(string filename) {
@ -22,6 +22,7 @@ void ScaLBL_Multiphys_Controller::ReadParams(string filename) {
visualization_interval = 10000;
tolerance = 1.0e-6;
time_conv_max = 0.0;
time_conv_MainLoop = 0.0;
// load input parameters
if (study_db->keyExists("timestepMax")) {

1
models/MultiPhysController.h
View File

@ -40,6 +40,7 @@ public:
int visualization_interval;
double tolerance;
double time_conv_max;
double time_conv_MainLoop;
//double SchmidtNum;//Schmidt number = kinematic_viscosity/mass_diffusivity
int rank, nprocs;

18
models/PoissonSolver.cpp
View File

@ -523,7 +523,7 @@ void ScaLBL_Poisson::Initialize(double time_conv_from_Study){
//}
}
void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study){
//.......create and start timer............
//double starttime,stoptime,cputime;
@ -537,13 +537,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
// *************ODD TIMESTEP*************//
timestep++;
SolveElectricPotentialAAodd(timestep_from_Study);//update electric potential
SolvePoissonAAodd(ChargeDensity);//perform collision
SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier();
// *************EVEN TIMESTEP*************//
timestep++;
SolveElectricPotentialAAeven(timestep_from_Study);//update electric potential
SolvePoissonAAeven(ChargeDensity);//perform collision
SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
@ -724,9 +724,9 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){
ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np);
}
void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np);
void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC){
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np);
//TODO: perhaps add another ScaLBL_Comm routine to update Psi values on solid boundary nodes.
//something like:
//ScaLBL_Comm->SolidDirichletBoundaryUpdates(Psi, Psi_BCLabel, timestep);
@ -739,9 +739,9 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){
//}
}
void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity){
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np);
void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC){
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->SolidDirichletAndNeumannD3Q7(fq, Psi, Psi_BCLabel);
//if (BoundaryConditionSolid==1){
// ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi);

6
models/PoissonSolver.h
View File

@ -34,7 +34,7 @@ public:
void ReadInput();
void Create();
void Initialize(double time_conv_from_Study);
void Run(double *ChargeDensity, int timestep_from_Study);
void Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study);
void getElectricPotential(DoubleArray &ReturnValues);
void getElectricPotential_debug(int timestep);
void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y,
@ -111,8 +111,8 @@ private:
void SolveElectricPotentialAAodd(int timestep_from_Study);
void SolveElectricPotentialAAeven(int timestep_from_Study);
//void SolveElectricField();
void SolvePoissonAAodd(double *ChargeDensity);
void SolvePoissonAAeven(double *ChargeDensity);
void SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC);
void SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC);
void getConvergenceLog(int timestep,double error);
double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int time_step);

14
models/StokesModel.cpp
View File

@ -593,7 +593,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
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,
rlx_setB, Fx, Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
@ -610,8 +610,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
}
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);
Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
0, ScaLBL_Comm->LastExterior(), Np);
if (UseSlippingVelBC == true) {
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(
@ -626,8 +626,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
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, ScaLBL_Comm->FirstInterior(),
ScaLBL_Comm->LastInterior(), Np);
Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
if (BoundaryCondition == 3) {
@ -643,8 +643,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
}
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);
Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
0, ScaLBL_Comm->LastExterior(), Np);
if (UseSlippingVelBC == true) {
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(
fq, ZetaPotentialSolid, ElectricField, SolidGrad, epsilon_LB,

2
sample_scripts/configure_arden
View File

@ -15,6 +15,6 @@ cmake -D CMAKE_C_COMPILER:PATH=/opt/arden/openmpi/3.1.2/bin/mpicc \
-D HDF5_DIRECTORY="/opt/arden/hdf5/1.8.12" \
-D USE_CUDA=0 \
-D USE_TIMER=0 \
~/Programs/LBPM
~/Programs/LBPM-WIA
# -D HDF5_LIB="/opt/arden/hdf5/1.8.12/lib/libhdf5.a"\

2
tests/TestNernstPlanck.cpp
View File

@ -74,7 +74,7 @@ int main(int argc, char **argv)
while (timestep < Study.timestepMax && error > Study.tolerance){
timestep++;
PoissonSolver.Run(IonModel.ChargeDensity,0);//solve Poisson equtaion to get steady-state electrical potental
PoissonSolver.Run(IonModel.ChargeDensity,false,0);//solve Poisson equtaion to get steady-state electrical potental
IonModel.Run(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField); //solve for ion transport and electric potential
timestep++;//AA operations

2
tests/TestPNP_Stokes.cpp
View File

@ -94,7 +94,7 @@ int main(int argc, char **argv)
while (timestep < Study.timestepMax && error > Study.tolerance){
timestep++;
PoissonSolver.Run(IonModel.ChargeDensity,0);//solve Poisson equtaion to get steady-state electrical potental
PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,0);//solve Poisson equtaion to get steady-state electrical potental
StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential

4
tests/TestPoissonSolver.cpp
View File

@ -69,7 +69,7 @@ int main(int argc, char **argv)
int timeSave = int(PoissonSolver.TestPeriodicSaveInterval/PoissonSolver.TestPeriodicTimeConv);
while (timestep<timeMax){
timestep++;
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,timestep);
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,false,timestep);
if (timestep%timeSave==0){
if (rank==0) printf(" Time = %.3g[s]; saving electric potential and field\n",timestep*PoissonSolver.TestPeriodicTimeConv);
PoissonSolver.getElectricPotential_debug(timestep);
@ -78,7 +78,7 @@ int main(int argc, char **argv)
}
}
else {
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,1);
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,false,1);
PoissonSolver.getElectricPotential_debug(1);
PoissonSolver.getElectricField_debug(1);
}

19
tests/lbpm_electrokinetic_SingleFluid_simulator.cpp
View File

@ -79,26 +79,35 @@ int main(int argc, char **argv)
IonModel.timestepMax = Study.getIonNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
IonModel.Initialize();
// Get maximal time converting factor based on Sotkes and Ion solvers
Study.getTimeConvMax_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
//Study.getTimeConvMax_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
// Get time conversion factor for the main iteration loop in electrokinetic single fluid simulator
Study.time_conv_MainLoop = StokesModel.timestepMax*StokesModel.time_conv;
// Initialize LB-Poisson model
PoissonSolver.ReadParams(filename);
PoissonSolver.SetDomain();
PoissonSolver.ReadInput();
PoissonSolver.Create();
PoissonSolver.Initialize(Study.time_conv_max);
PoissonSolver.Initialize(Study.time_conv_MainLoop);
printf("********************************************************\n");
printf("Key Summary of LBPM electrokinetic single-fluid solver \n");
printf(" 1. Max LB Timestep: %i [lt]\n", Study.timestepMax);
printf(" 2. Time conversion factor per LB Timestep: %.6g [sec/lt]\n",Study.time_conv_MainLoop);
printf(" 3. Max Physical Time: %.6g [sec]\n",Study.timestepMax*Study.time_conv_MainLoop);
printf("********************************************************\n");
int timestep=0;
while (timestep < Study.timestepMax){
timestep++;
PoissonSolver.Run(IonModel.ChargeDensity,timestep);//solve Poisson equtaion to get steady-state electrical potental
PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
timestep++;//AA operations
if (timestep%Study.analysis_interval==0){
Analysis.Basic(IonModel,PoissonSolver,StokesModel,timestep);
@ -116,7 +125,7 @@ int main(int argc, char **argv)
if (rank==0) printf("Save simulation raw data at maximum timestep\n");
Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep);
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
if (rank==0) printf("Maximum LB timestep = %i is reached and the simulation is completed\n",Study.timestepMax);
if (rank==0) printf("*************************************************************\n");
PROFILE_STOP("Main");

3
tests/lbpm_morphdrain_pp.cpp
View File

@ -53,6 +53,7 @@ int main(int argc, char **argv)
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
SW = domain_db->getScalar<double>("Sw");
auto READFILE = domain_db->getScalar<std::string>( "Filename" );
auto MORPH_RADIUS = domain_db->getWithDefault<double>( "MorphRadius", 100000.0);
// Generate the NWP configuration
//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
@ -122,7 +123,7 @@ int main(int argc, char **argv)
comm.barrier();
// Run the morphological opening
MorphDrain(SignDist, id, Dm, SW);
MorphDrain(SignDist, id, Dm, SW, MORPH_RADIUS);
// calculate distance to non-wetting fluid
if (domain_db->keyExists( "HistoryLabels" )){