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

save the work;untested

Browse Source
This commit is contained in:
Rex Zhe Li
2020-08-10 12:03:28 -04:00
parent dff4e3d536
commit 28988ef6ba
5 changed files with 395 additions and 127 deletions
Download Patch File Download Diff File Expand all files Collapse all files

316
cpu/Ion.cpp
View File

@@ -1,12 +1,177 @@
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Velocity, 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;
for (int n=start; n<finish; n++){ extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *dist, double *Den, int start, int finish, int Np){
int n,nread;
double fq,Ci;
for (n=start; n<finish; n++){
// q=0 // q=0
fq = dist[n];
Ci = fq;
// q=1
nread = neighborList[n];
fq = dist[nread];
Ci += fq;
// q=2
nread = neighborList[n+Np];
fq = dist[nread];
Ci += fq;
// q=3
nread = neighborList[n+2*Np];
fq = dist[nread];
Ci += fq;
// q=4
nread = neighborList[n+3*Np];
fq = dist[nread];
Ci += fq;
// q=5
nread = neighborList[n+4*Np];
fq = dist[nread];
Ci += fq;
// q=6
nread = neighborList[n+5*Np];
fq = dist[nread];
Ci += fq;
Den[n]=Ci;
}
}
extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, int start, int finish, int Np){
int n;
double fq,Ci;
for (n=start; n<finish; n++){
// q=0
fq = dist[n];
Ci = fq;
// q=1
fq = dist[2*Np+n];
Ci += fq;
// q=2
fq = dist[1*Np+n];
Ci += fq;
// q=3
fq = dist[4*Np+n];
Ci += fq;
// q=4
fq = dist[3*Np+n];
Ci += fq;
// q=5
fq = dist[6*Np+n];
Ci += fq;
// q=6
fq = dist[5*Np+n];
Ci += fq;
Den[n]=Ci;
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, double zi, double rlx, double deltaT, 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 f0,f1,f2,f3,f4,f5,f6;
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];
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*Ci;
// q = 1
dist[nr2] = f1*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(ux+uEPx));
// q=2
dist[nr1] = f2*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(ux+uEPx));
// q = 3
dist[nr4] = f3*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(uy+uEPy));
// q = 4
dist[nr3] = f4*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(uy+uEPy));
// q = 5
dist[nr6] = f5*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(uz+uEPz));
// q = 6
dist[nr5] = f6*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(uz+uEPz));
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, double zi, double rlx, double deltaT, 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 f0,f1,f2,f3,f4,f5,f6;
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]; f0 = dist[n];
f1 = dist[2*Np+n]; f1 = dist[2*Np+n];
f2 = dist[1*Np+n]; f2 = dist[1*Np+n];
@@ -14,110 +179,67 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Velocity, int start
f4 = dist[3*Np+n]; f4 = dist[3*Np+n];
f5 = dist[6*Np+n]; f5 = dist[6*Np+n];
f6 = dist[5*Np+n]; f6 = dist[5*Np+n];
rho = f0+f2+f1+f4+f3+f6;
ux = Velocity[n];
uy = Velocity[n+Np];
uz = Velocity[n+2*Np];
uu = 1.5*(ux*ux+uy*uy+uz*uz);
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*(1.0-uu);
// q = 1
dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.05555555555555555*(rho + 3.0*ux + 4.5*ux*ux - uu) + 0.16666666*Fx;
// q=2
dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.05555555555555555*(rho - 3.0*ux + 4.5*ux*ux - uu)- 0.16666666*Fx;
// q = 3
dist[3*Np+n] = f3*(1.0-rlx) +
rlx*0.05555555555555555*(rho + 3.0*uy + 4.5*uy*uy - uu) + 0.16666666*Fy;
// q = 4
dist[4*Np+n] = f4*(1.0-rlx) +
rlx*0.05555555555555555*(rho - 3.0*uy + 4.5*uy*uy - uu)- 0.16666666*Fy;
// q = 5
dist[5*Np+n] = f5*(1.0-rlx) +
rlx*0.05555555555555555*(rho + 3.0*uz + 4.5*uz*uz - uu) + 0.16666666*Fz;
// q = 6
dist[6*Np+n] = f6*(1.0-rlx) +
rlx*0.05555555555555555*(rho - 3.0*uz + 4.5*uz*uz - uu) - 0.16666666*Fz;
//........................................................................
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Velocity, 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 // q=0
f0 = dist[n]; dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*Ci;
// q=1
nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
f1 = dist[nr1]; // reading the f1 data into register fq
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];
rho = f0+f2+f1+f4+f3+f6;
ux = Velocity[n];
uy = Velocity[n+Np];
uz = Velocity[n+2*Np];
uu = 1.5*(ux*ux+uy*uy+uz*uz);
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*(1.0-uu);
// q = 1 // q = 1
dist[nr2] = f1*(1.0-rlx) + rlx*0.05555555555555555*(rho + 3.0*ux + 4.5*ux*ux - uu) + 0.16666666*Fx; dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(ux+uEPx));
// q=2 // q=2
dist[nr1] = f2*(1.0-rlx) + rlx*0.05555555555555555*(rho - 3.0*ux + 4.5*ux*ux - uu)- 0.16666666*Fx; dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(ux+uEPx));
// q = 3 // q = 3
dist[nr4] = f3*(1.0-rlx) + dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(uy+uEPy));
rlx*0.05555555555555555*(rho + 3.0*uy + 4.5*uy*uy - uu) + 0.16666666*Fy;
// q = 4 // q = 4
dist[nr3] = f4*(1.0-rlx) + dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(uy+uEPy));
rlx*0.05555555555555555*(rho - 3.0*uy + 4.5*uy*uy - uu)- 0.16666666*Fy;
// q = 5 // q = 5
dist[nr6] = f5*(1.0-rlx) + dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(uz+uEPz));
rlx*0.05555555555555555*(rho + 3.0*uz + 4.5*uz*uz - uu) + 0.16666666*Fz;
// q = 6 // q = 6
dist[nr5] = f6*(1.0-rlx) + dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(uz+uEPz));
rlx*0.05555555555555555*(rho - 3.0*uz + 4.5*uz*uz - uu) - 0.16666666*Fz;
} }
} }
extern "C" void ScaLBL_D3Q7_Poisson_Init(double *dist, double *Den, double DenInit, int Np)
{
int n;
for (n=0; n<Np; n++){
dist[0*Np+n] = 0.3333333333333333*DenInit;
dist[1*Np+n] = 0.1111111111111111*DenInit;
dist[2*Np+n] = 0.1111111111111111*DenInit;
dist[3*Np+n] = 0.1111111111111111*DenInit;
dist[4*Np+n] = 0.1111111111111111*DenInit;
dist[5*Np+n] = 0.1111111111111111*DenInit;
dist[6*Np+n] = 0.1111111111111111*DenInit;
Den[n] = DenInit;
}
}
extern "C" void ScaLBL_D3Q7_IonChargeDensity(double *Den, double *ChargeDensity, vector<double>& IonValence, int number_ion_species, int start, int finish, int Np){
int n;
int ic=number_ion_species;
double Ci;//ion concentration of species i
double CD;//charge density
double F = 96485.0;//Faraday's constant; unit[C/mol]; F=e*Na, where Na is the Avogadro constant
for (n=start; n<finish; n++){
Ci = Den[n];
CD = F*IonValence[0]*Ci;
ChargeDensity[n] = CD;
}
ic = ic - 1;
while (ic>0){
for (n=start; n<finish; n++){
Ci = Den[n+ic*Np];
CD = F*IonValence[ic]*Ci;
ChargeDensity[n] += CD;
}
ic--;
}
}

4
cpu/Poisson.cpp
View File

@@ -8,7 +8,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, doubl
double f0,f1,f2,f3,f4,f5,f6; double f0,f1,f2,f3,f4,f5,f6;
int nr1,nr2,nr3,nr4,nr5,nr6; int nr1,nr2,nr3,nr4,nr5,nr6;
for (int n=start; n<finish; n++){ for (n=start; n<finish; n++){
//Load data //Load data
rho_e = Den_charge[n]; rho_e = Den_charge[n];
@@ -80,7 +80,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, dou
double rho_e;//local charge density double rho_e;//local charge density
double f0,f1,f2,f3,f4,f5,f6; double f0,f1,f2,f3,f4,f5,f6;
for (int n=start; n<finish; n++){ for (n=start; n<finish; n++){
//Load data //Load data
rho_e = Den_charge[n]; rho_e = Den_charge[n];

186
models/IonModel.cpp
View File

@@ -17,21 +17,103 @@ ScaLBL_IonModel::~ScaLBL_IonModel(){
} }
void ScaLBL_IonModel::ReadParams(string filename){ void ScaLBL_IonModel::ReadParams(string filename){
// read the input database
//fundamental constant
kb = 1.38e-23;//Boltzmann constant;unit [J/K]
electron_charge = 1.6e-19;//electron charge;unit [C]
// read the input database
db = std::make_shared<Database>( filename ); db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" ); domain_db = db->getDatabase( "Domain" );
ion_db = db->getDatabase( "Ions" ); ion_db = db->getDatabase( "Ions" );
tau = 1.0; // Default model parameters
T = 300.0;//temperature; unit [K]
Vt = kb*T/electron_charge;//thermal voltage; unit [V]
k2_inv = 4.5;//the inverse of 2nd-rank moment of D3Q7 lattice
h = 1.0;//resolution; unit: um/lu
timestepMax = 100000; timestepMax = 100000;
tolerance = 1.0e-8; tolerance = 1.0e-8;
// Model parameters
number_ion_species = 1; number_ion_species = 1;
IonDiffusivity.push_back(1.0e-9);//User input unit [m^2/sec]
//TODO needs to scale the unit of diffusivity!
IonValence.push_back(1);
IonConcentration.push_back(1.0e-3);//unit [mol/m^3]
// TODO rescale ion concentration unit
deltaT.push_back(1.0);
tau.push_back(0.5+k2_inv*deltaT[0]*IonDiffusivisty[0]);
// LB-Ion Model parameters
if (ion_db->keyExists( "timestepMax" )){
timestepMax = ion_db->getScalar<int>( "timestepMax" );
}
if (ion_db->keyExists( "analysis_interval" )){
analysis_interval = ion_db->getScalar<int>( "analysis_interval" );
}
if (ion_db->keyExists( "tolerance" )){
tolerance = ion_db->getScalar<double>( "tolerance" );
}
if (ion_db->keyExists( "temperature" )){
T = ion_db->getScalar<int>( "temperature" );
}
if (ion_db->keyExists( "epsilonR" )){
epsilonR = ion_db->getScalar<double>( "epsilonR" );
}
if (ion_db->keyExists( "number_ion_species" )){ if (ion_db->keyExists( "number_ion_species" )){
number_ion_species = ion_db->getScalar<int>( "number_ion_species" ); number_ion_species = ion_db->getScalar<int>( "number_ion_species" );
} }
//read ion related list
if (ion_db->keyExists( "deltaT" )){
deltaT.clear();
deltaT = ion_db->getVector<double>( "deltaT" );
if (deltaT.size()!=number_ion_species){
ERROR("Error: number_ion_species and deltaT must be the same length! \n");
}
}
//NOTE: Ion diffusivity has unit: [m^2/sec]
if (ion_db->keyExists("IonDiffusivityList")){
IonDiffusivity.clear();
IonDiffusivity = ion_db->getVector<double>( "IonDiffusivityList" );
if (IonDiffusivity.size()!=number_ion_species){
ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n");
}
}
//read ion algebric valence list
if (ion_db->keyExists("IonValenceList")){
IonValence.clear();
IonValence = ion_db->getVector<int>( "IonValenceList" );
if (IonValence.size()!=number_ion_species){
ERROR("Error: number_ion_species and IonValenceList must be the same length! \n");
}
}
//read initial ion concentration list; unit [mol/m^3]
if (ion_db->keyExists("IonConcentrationList")){
IonConcentration.clear();
IonConcentration = ion_db->getVector<double>( "IonConcentrationList" );
if (IonConcentration.size()!=number_ion_species){
ERROR("Error: number_ion_species and IonConcentrationList must be the same length! \n");
}
}
// Read domain parameters
if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
h = domain_db->getScalar<double>( "voxel_length" );
}
if (domain_db->keyExists( "BC" )){
BoundaryCondition = domain_db->getScalar<int>( "BC" );
}
//Re-calcualte model parameters if user updates input
//TODO ion diffusivity needs rescale unit to LB unit
//TODO rescale ion initial concentration unit to LB unit
if (deltaT.size()>1){
tau.clear();
for (int i=0;i<IonDiffusivity.size();i++){
tau.push_back(0.5+k2_inv*deltaT*IonDiffusivity[i]);
}
}
} }
void ScaLBL_IonModel::SetDomain(){ void ScaLBL_IonModel::SetDomain(){
Dm = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // full domain for analysis Dm = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // full domain for analysis
Mask = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // mask domain removes immobile phases Mask = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // mask domain removes immobile phases
@@ -164,42 +246,104 @@ void ScaLBL_IonModel::Initialize(){
/* /*
* This function initializes model * This function initializes model
*/ */
if (rank==0) printf ("Initializing distributions \n"); if (rank==0) printf ("Initializing D3Q7 distributions for ion transport\n");
ScaLBL_D3Q19_Init(fq, Np); for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init(&fq[ic*Np*7],&Ci[ic*Np],IonConcentration[ic],Np);
}
} }
void ScaLBL_IonModel::Run(double *Velocity){ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
double rlx = 1.0/tau;
//LB-related parameter
vector<double> rlx(tau.begin(),tau.end());
for (double item : rlx){
item = 1.0/item;
}
//.......create and start timer............ //.......create and start timer............
double starttime,stoptime,cputime; double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier(); MPI_Barrier(comm); ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
starttime = MPI_Wtime(); starttime = MPI_Wtime();
if (rank==0) printf("Beginning AA timesteps, timestepMax = %i \n", timestepMax);
if (rank==0) printf("********************************************************\n"); if (rank==0) printf("***************************************************\n");
if (rank==0) printf("LB-Ion Transport: timestepMax = %i\n", timestepMax);
if (rank==0) printf("***************************************************\n");
timestep=0; timestep=0;
double error = 1.0; while (timestep < timestepMax) {
double flow_rate_previous = 0.0;
while (timestep < timestepMax && error > tolerance) {
//************************************************************************/ //************************************************************************/
// *************ODD TIMESTEP*************//
timestep++; timestep++;
for (int ic=0; ic<number_ion_species; ic++) //Update ion concentration and charge density
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FROM NORMAL ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_Ion(NeighborList, fq, Velocity, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz); }
for (int ic=0; ic<number_ion_species; ic++) for (int ic=0; ic<number_ion_species; ic++){
//TODO should this loop be merged with Send & Recv ?
//Sum up distribution to get ion concentration
ScaLBL_D3Q7_AAodd_IonConcentration(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
ScaLBL_D3Q7_IonChargeDensity(Ci, ChargeDensity, IonValence, number_ion_species, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE
}
for (int ic=0; ic<number_ion_species; ic++){
//TODO should this loop be merged with Send & Recv ?
//Sum up distribution to get ion concentration
ScaLBL_D3Q7_AAodd_IonConcentration(NeighborList, &fq[ic*Np*7],&Ci[ic*Np], 0, ScaLBL_Comm->LastExterior(), Np);
}
ScaLBL_D3Q7_IonChargeDensity(Ci, ChargeDensity, number_ion_species, 0, ScaLBL_Comm->LastExterior(), Np);
//LB-Ion collison
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],deltaT[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
// Set boundary conditions // Set boundary conditions
/* ... */ /* ... */
ScaLBL_D3Q7_AAodd_Ion(NeighborList, fq, Velocity, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],deltaT[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm); ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// *************EVEN TIMESTEP*************//
timestep++; timestep++;
for (int ic=0; ic<number_ion_species; ic++) //Update ion concentration and charge density
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FORM NORMAL ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FORM NORMAL
ScaLBL_D3Q7_AAeven_Ion(fq, Velocity, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz); }
for (int ic=0; ic<number_ion_species; ic++) for (int ic=0; ic<number_ion_species; ic++){
//TODO should this loop be merged with Send & Recv ?
//Sum up distribution to get ion concentration
ScaLBL_D3Q7_AAeven_IonConcentration(&fq[ic*Np*7],&Ci[ic*Np],ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
ScaLBL_D3Q7_IonChargeDensity(Ci, ChargeDensity, IonValence, number_ion_species, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE
}
for (int ic=0; ic<number_ion_species; ic++){
//TODO should this loop be merged with Send & Recv ?
//Sum up distribution to get ion concentration
ScaLBL_D3Q7_AAeven_IonConcentration(&fq[ic*Np*7],&Ci[ic*Np], 0, ScaLBL_Comm->LastExterior(), Np);
}
ScaLBL_D3Q7_IonChargeDensity(Ci, ChargeDensity, number_ion_species, 0, ScaLBL_Comm->LastExterior(), Np);
//LB-Ion collison
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],deltaT[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
// Set boundary conditions // Set boundary conditions
/* ... */ /* ... */
ScaLBL_D3Q7_AAeven_Ion(fq, Velocity, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],deltaT[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm); ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/ //************************************************************************/
} }

10
models/IonModel.h
View File

@@ -8,6 +8,7 @@
#include <exception> #include <exception>
#include <stdexcept> #include <stdexcept>
#include <fstream> #include <fstream>
#include <vector>
#include "common/ScaLBL.h" #include "common/ScaLBL.h"
#include "common/Communication.h" #include "common/Communication.h"
@@ -33,12 +34,13 @@ public:
bool Restart,pBC; bool Restart,pBC;
int timestep,timestepMax; int timestep,timestepMax;
int BoundaryCondition; int BoundaryCondition;
double tau,mu;
double Fx,Fy,Fz,flux;
double din,dout;
double tolerance;
int number_ion_species; int number_ion_species;
vector<double> IonDiffusivity;//User input unit [m^2/sec]
vector<int> IonValence;
vector<double> IonConcentration;//unit [mol/m^3]
vector<double> deltaT;
vector<double> tau;
int Nx,Ny,Nz,N,Np; int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs; int rank,nprocx,nprocy,nprocz,nprocs;

6
tests/lbpm_electrokinetic_dfh_simulator.cpp
View File

@@ -59,10 +59,10 @@ int main(int argc, char **argv)
DFHModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables DFHModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables
DFHModel.Initialize(); // initializing the model will set initial conditions for variables DFHModel.Initialize(); // initializing the model will set initial conditions for variables
DFHModel.Run(); // Solve the N-S equations to get velocity DFHModel.Run(IonModel.ChargeDensity); // Solve the N-S equations to get velocity
IonModel.Run(DFHModel.Velocity); //solve for ion transport and electric potential
//IonModel.Run(DFHModel.Velocity, DFHModel.Phi); //solve for ion transport and electric potential with multiphase system
PoissonSolver.Run(IonModel.ChargeDensity); PoissonSolver.Run(IonModel.ChargeDensity);
IonModel.Run(DFHModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
//IonModel.Run(DFHModel.Velocity, DFHModel.Phi,PoissonSolver.ElectricField);
DFHModel.WriteDebug(); DFHModel.WriteDebug();