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

correct ion LB collison

Browse Source
This commit is contained in:
Zhe Rex Li 2022-05-03 15:00:00 +10:00
parent 47db000ba6
commit 6959a02f7b
3 changed files with 220 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
common/ScaLBL.h
View File

@ -291,6 +291,11 @@ extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *di
extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion_v0(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Ion_v0(double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective, double *FluxElectrical, double *Velocity, double *ElectricField,
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);

213
cpu/Ion.cpp
View File

@ -242,7 +242,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den,
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
extern "C" void ScaLBL_D3Q7_AAodd_Ion_v0(int *neighborList, double *dist,
double *Den, double *FluxDiffusive,
double *FluxAdvective,
double *FluxElectrical, double *Velocity,
@ -356,7 +356,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Ion(
extern "C" void ScaLBL_D3Q7_AAeven_Ion_v0(
double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective,
double *FluxElectrical, double *Velocity, double *ElectricField, double Di,
int zi, double rlx, double Vt, int start, int finish, int Np) {
@ -451,6 +451,215 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
double *Den, double *FluxDiffusive,
double *FluxAdvective,
double *FluxElectrical, double *Velocity,
double *ElectricField, double Di, int zi,
double rlx, double Vt, int start,
int finish, int Np) {
int n;
double Ci;
double ux, uy, uz;
double uEPx, uEPy, uEPz; //electrochemical induced velocity
double Ex, Ey, Ez; //electrical field
double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z;
double f0, f1, f2, f3, f4, f5, f6;
//double X,Y,Z,factor_x, factor_y, factor_z;
int nr1, nr2, nr3, nr4, nr5, nr6;
for (n = start; n < finish; n++) {
//Load data
//Ci = Den[n];
Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
ux = Velocity[n + 0 * Np];
uy = Velocity[n + 1 * Np];
uz = Velocity[n + 2 * Np];
uEPx = zi * Di / Vt * Ex;
uEPy = zi * Di / Vt * Ey;
uEPz = zi * Di / Vt * Ez;
// q=0
f0 = dist[n];
// q=1
nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
f1 = dist[nr1]; // reading the f1 data into register fq
// q=2
nr2 = neighborList[n + Np]; // neighbor 1 ( < 10Np => even part of dist)
f2 = dist[nr2]; // reading the f2 data into register fq
// q=3
nr3 = neighborList[n + 2 * Np]; // neighbor 4
f3 = dist[nr3];
// q=4
nr4 = neighborList[n + 3 * Np]; // neighbor 3
f4 = dist[nr4];
// q=5
nr5 = neighborList[n + 4 * Np];
f5 = dist[nr5];
// q=6
nr6 = neighborList[n + 5 * Np];
f6 = dist[nr6];
// compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
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_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
FluxDiffusive[n + 0 * Np] = flux_diffusive_x;
FluxDiffusive[n + 1 * Np] = flux_diffusive_y;
FluxDiffusive[n + 2 * Np] = flux_diffusive_z;
FluxAdvective[n + 0 * Np] = ux * Ci;
FluxAdvective[n + 1 * Np] = uy * Ci;
FluxAdvective[n + 2 * Np] = uz * Ci;
FluxElectrical[n + 0 * Np] = uEPx * Ci;
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci;
/* use logistic function to prevent negative distributions*/
//X = 4.0 * (ux + uEPx);
//Y = 4.0 * (uy + uEPy);
//Z = 4.0 * (uz + uEPz);
//factor_x = X / sqrt(1 + X*X);
//factor_y = Y / sqrt(1 + Y*Y);
//factor_z = Z / sqrt(1 + Z*Z);
// q=0
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1
dist[nr2] =
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
dist[nr1] =
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
dist[nr4] =
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
dist[nr3] =
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
dist[nr6] =
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
dist[nr5] =
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);
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Ion(
double *dist, double *Den, double *FluxDiffusive, double *FluxAdvective,
double *FluxElectrical, double *Velocity, double *ElectricField, double Di,
int zi, double rlx, double Vt, int start, int finish, int Np) {
int n;
double Ci;
double ux, uy, uz;
double uEPx, uEPy, uEPz; //electrochemical induced velocity
double Ex, Ey, Ez; //electrical field
double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z;
double f0, f1, f2, f3, f4, f5, f6;
//double X,Y,Z, factor_x, factor_y, factor_z;
for (n = start; n < finish; n++) {
//Load data
//Ci = Den[n];
Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
ux = Velocity[n + 0 * Np];
uy = Velocity[n + 1 * Np];
uz = Velocity[n + 2 * Np];
uEPx = zi * Di / Vt * Ex;
uEPy = zi * Di / Vt * Ey;
uEPz = zi * Di / Vt * Ez;
f0 = dist[n];
f1 = dist[2 * Np + n];
f2 = dist[1 * Np + n];
f3 = dist[4 * Np + n];
f4 = dist[3 * Np + n];
f5 = dist[6 * Np + n];
f6 = dist[5 * Np + n];
// compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
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_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
FluxDiffusive[n + 0 * Np] = flux_diffusive_x;
FluxDiffusive[n + 1 * Np] = flux_diffusive_y;
FluxDiffusive[n + 2 * Np] = flux_diffusive_z;
FluxAdvective[n + 0 * Np] = ux * Ci;
FluxAdvective[n + 1 * Np] = uy * Ci;
FluxAdvective[n + 2 * Np] = uz * Ci;
FluxElectrical[n + 0 * Np] = uEPx * Ci;
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci;
/* use logistic function to prevent negative distributions*/
//X = 4.0 * (ux + uEPx);
//Y = 4.0 * (uy + uEPy);
//Z = 4.0 * (uz + uEPz);
//factor_x = X / sqrt(1 + X*X);
//factor_y = Y / sqrt(1 + Y*Y);
//factor_z = Z / sqrt(1 + Z*Z);
// q=0
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1
dist[1 * Np + n] =
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
dist[2 * Np + n] =
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
dist[3 * Np + n] =
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
dist[4 * Np + n] =
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
dist[5 * Np + n] =
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
dist[6 * Np + n] =
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);
}
}
extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit,
int Np) {
int n;

8
models/IonModel.cpp
View File

@ -1296,13 +1296,13 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField) {
ScaLBL_Comm->LastExterior(), Np);
//LB-Ion collison
ScaLBL_D3Q7_AAodd_Ion(
ScaLBL_D3Q7_AAodd_Ion_v0(
NeighborList, &fq[ic * Np * 7], &Ci[ic * Np],
&FluxDiffusive[3 * ic * Np], &FluxAdvective[3 * ic * Np],
&FluxElectrical[3 * ic * Np], Velocity, ElectricField,
IonDiffusivity[ic], IonValence[ic], rlx[ic], Vt,
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_Ion(
ScaLBL_D3Q7_AAodd_Ion_v0(
NeighborList, &fq[ic * Np * 7], &Ci[ic * Np],
&FluxDiffusive[3 * ic * Np], &FluxAdvective[3 * ic * Np],
&FluxElectrical[3 * ic * Np], Velocity, ElectricField,
@ -1380,13 +1380,13 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField) {
Np);
//LB-Ion collison
ScaLBL_D3Q7_AAeven_Ion(
ScaLBL_D3Q7_AAeven_Ion_v0(
&fq[ic * Np * 7], &Ci[ic * Np], &FluxDiffusive[3 * ic * Np],
&FluxAdvective[3 * ic * Np], &FluxElectrical[3 * ic * Np],
Velocity, ElectricField, IonDiffusivity[ic], IonValence[ic],
rlx[ic], Vt, ScaLBL_Comm->FirstInterior(),
ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_Ion(
ScaLBL_D3Q7_AAeven_Ion_v0(
&fq[ic * Np * 7], &Ci[ic * Np], &FluxDiffusive[3 * ic * Np],
&FluxAdvective[3 * ic * Np], &FluxElectrical[3 * ic * Np],
Velocity, ElectricField, IonDiffusivity[ic], IonValence[ic],