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

make further corrections

Browse Source
This commit is contained in:
Zhe Li
2022-04-28 23:59:13 +10:00
parent bdd0efd36e
commit 025987e53b
2 changed files with 45 additions and 43 deletions
Download Patch File Download Diff File Expand all files Collapse all files

86
cpu/Ion.cpp
View File

@@ -256,12 +256,13 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
double Ex, Ey, Ez; //electrical field double Ex, Ey, Ez; //electrical field
double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z; double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z;
double f0, f1, f2, f3, f4, f5, f6; double f0, f1, f2, f3, f4, f5, f6;
double X,Y,Z,factor_x, factor_y, factor_z; //double X,Y,Z,factor_x, factor_y, factor_z;
int nr1, nr2, nr3, nr4, nr5, nr6; int nr1, nr2, nr3, nr4, nr5, nr6;
for (n = start; n < finish; n++) { for (n = start; n < finish; n++) {
//Load data //Load data
Ci = Den[n];
Ex = ElectricField[n + 0 * Np]; Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np]; Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np]; Ez = ElectricField[n + 2 * Np];
@@ -294,7 +295,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
f6 = dist[nr6]; f6 = dist[nr6];
// compute diffusive flux // compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6; //Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci); flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci); flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci); flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@@ -308,48 +309,49 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
FluxElectrical[n + 1 * Np] = uEPy * Ci; FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci; FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci; //Den[n] = Ci;
/* use logistic function to prevent negative distributions*/ /* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx); //X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy); //Y = 4.0 * (uy + uEPy);
Z = 4.0 * (uz + uEPz); //Z = 4.0 * (uz + uEPz);
factor_x = X / sqrt(1 + X*X); //factor_x = X / sqrt(1 + X*X);
factor_y = Y / sqrt(1 + Y*Y); //factor_y = Y / sqrt(1 + Y*Y);
factor_z = Z / sqrt(1 + Z*Z); //factor_z = Z / sqrt(1 + Z*Z);
// q=0 // q=0
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci; dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1 // q = 1
dist[nr2] = dist[nr2] =
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x); f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
//f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx)); // f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
// q=2 // q=2
dist[nr1] = dist[nr1] =
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x); f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
//f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx)); // f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
// q = 3 // q = 3
dist[nr4] = dist[nr4] =
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y ); f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
//f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy)); // f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y );
// q = 4 // q = 4
dist[nr3] = dist[nr3] =
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y); f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
//f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy)); // f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
// q = 5 // q = 5
dist[nr6] = dist[nr6] =
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z); f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
//f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz)); // f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
// q = 6 // q = 6
dist[nr5] = dist[nr5] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z); f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
// f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
} }
} }
@@ -365,12 +367,12 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
double Ex, Ey, Ez; //electrical field double Ex, Ey, Ez; //electrical field
double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z; double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z;
double f0, f1, f2, f3, f4, f5, f6; double f0, f1, f2, f3, f4, f5, f6;
double X,Y,Z, factor_x, factor_y, factor_z; //double X,Y,Z, factor_x, factor_y, factor_z;
for (n = start; n < finish; n++) { for (n = start; n < finish; n++) {
//Load data //Load data
//Ci = Den[n]; Ci = Den[n];
Ex = ElectricField[n + 0 * Np]; Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np]; Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np]; Ez = ElectricField[n + 2 * Np];
@@ -390,7 +392,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
f6 = dist[5 * Np + n]; f6 = dist[5 * Np + n];
// compute diffusive flux // compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6; //Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci); flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci); flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci); flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@@ -404,48 +406,48 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
FluxElectrical[n + 1 * Np] = uEPy * Ci; FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci; FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci; //Den[n] = Ci;
/* use logistic function to prevent negative distributions*/ /* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx); //X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy); //Y = 4.0 * (uy + uEPy);
Z = 4.0 * (uz + uEPz); //Z = 4.0 * (uz + uEPz);
factor_x = X / sqrt(1 + X*X); //factor_x = X / sqrt(1 + X*X);
factor_y = Y / sqrt(1 + Y*Y); //factor_y = Y / sqrt(1 + Y*Y);
factor_z = Z / sqrt(1 + Z*Z); //factor_z = Z / sqrt(1 + Z*Z);
// q=0 // q=0
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci; dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1 // q = 1
dist[1 * Np + n] = dist[1 * Np + n] =
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x); f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
//f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx)); // f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
// q=2 // q=2
dist[2 * Np + n] = dist[2 * Np + n] =
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x); f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
//f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx)); // f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
// q = 3 // q = 3
dist[3 * Np + n] = dist[3 * Np + n] =
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y); f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
//f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy)); // f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y);
// q = 4 // q = 4
dist[4 * Np + n] = dist[4 * Np + n] =
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y); f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
//f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy)); // f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
// q = 5 // q = 5
dist[5 * Np + n] = dist[5 * Np + n] =
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z); f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
//f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz)); // f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
// q = 6 // q = 6
dist[6 * Np + n] = dist[6 * Np + n] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z); f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
//f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz)); // f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
} }
} }

2
tests/lbpm_nernst_planck_simulator.cpp
View File

@@ -90,7 +90,7 @@ int main(int argc, char **argv)
if (rank == 0){ if (rank == 0){
printf("********************************************************\n"); printf("********************************************************\n");
printf("Key Summary of LBPM electrokinetic single-fluid solver \n"); printf("Key Summary of LBPM electrokinetic single-fluid solver \n");
printf(" 1. Max LB Timestep: %i [lt]\n", Study.timestepMax); printf(" 1. Max LB Timestep: %li [lt]\n", Study.timestepMax);
printf(" 2. Time conversion factor per LB Timestep: %.6g [sec/lt]\n",Study.time_conv_MainLoop); 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(" 3. Max Physical Time: %.6g [sec]\n",Study.timestepMax*Study.time_conv_MainLoop);
printf("********************************************************\n"); printf("********************************************************\n");