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-14 14:23:22 -04:00
parent b4d50ee821
commit 3a162849a0
11 changed files with 1932 additions and 222 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
cpu/Ion.cpp
View File

@ -80,7 +80,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i
}
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){
double Di, double zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
@ -128,29 +128,28 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
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));
dist[nr2] = f1*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(ux+uEPx));
// q=2
dist[nr1] = f2*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(ux+uEPx));
dist[nr1] = f2*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(ux+uEPx));
// q = 3
dist[nr4] = f3*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(uy+uEPy));
dist[nr4] = f3*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(uy+uEPy));
// q = 4
dist[nr3] = f4*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(uy+uEPy));
dist[nr3] = f4*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(uy+uEPy));
// q = 5
dist[nr6] = f5*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(uz+uEPz));
dist[nr6] = f5*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(uz+uEPz));
// q = 6
dist[nr5] = f6*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(uz+uEPz));
dist[nr5] = f6*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(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){
double Di, double zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
@ -184,28 +183,28 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Veloci
dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*Ci;
// q = 1
dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(ux+uEPx));
dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(ux+uEPx));
// q=2
dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(ux+uEPx));
dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(ux+uEPx));
// q = 3
dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(uy+uEPy));
dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(uy+uEPy));
// q = 4
dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(uy+uEPy));
dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(uy+uEPy));
// q = 5
dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*deltaT*(uz+uEPz));
dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.1111111111111111*(1.0+4.5*(uz+uEPz));
// q = 6
dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*deltaT*(uz+uEPz));
dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.1111111111111111*(1.0-4.5*(uz+uEPz));
}
}
extern "C" void ScaLBL_D3Q7_Poisson_Init(double *dist, double *Den, double DenInit, int Np)
extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np)
{
int n;
for (n=0; n<Np; n++){
@ -220,7 +219,7 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(double *dist, double *Den, double DenIn
}
}
extern "C" void ScaLBL_D3Q7_IonChargeDensity(double *Den, double *ChargeDensity, vector<double>& IonValence, int number_ion_species, int start, int finish, int Np){
extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, vector<double>& IonValence, int number_ion_species, int start, int finish, int Np){
int n;
int ic=number_ion_species;

22
cpu/Poisson.cpp
View File

@ -1,5 +1,5 @@
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, double *Den_charge, double *Psi, double *ElectricField, double rlx, double epsilon_LB,double deltaT,
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,
int start, int finish, int Np){
int n;
double psi;//electric potential
@ -7,12 +7,13 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, doubl
double rho_e;//local charge density
double f0,f1,f2,f3,f4,f5,f6;
int nr1,nr2,nr3,nr4,nr5,nr6;
double rlx=1.0/tau;
for (n=start; n<finish; n++){
//Load data
rho_e = Den_charge[n];
rho_e = deltaT*rho_e/epsilon_LB;
rho_e = gamma*rho_e/epsilon_LB;
// q=0
f0 = dist[n];
@ -40,9 +41,9 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, doubl
f6 = dist[nr6];
psi = f0+f2+f1+f4+f3+f6+f5;
Ex = f1-f2;
Ey = f3-f4;
Ez = f5-f6;
Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
Ey = (f3-f4)*rlx*4.5;
Ez = (f5-f6)*rlx*4.5;
ElectricField[n+0*Np] = Ex;
ElectricField[n+1*Np] = Ey;
ElectricField[n+2*Np] = Ez;
@ -72,19 +73,20 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, doubl
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, double *Psi, double *ElectricField, double rlx, double epsilon_LB,double deltaT,
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,double gamma,
int start, int finish, int Np){
int n;
double psi;//electric potential
double Ex,Ey,Ez;//electrical field
double rho_e;//local charge density
double f0,f1,f2,f3,f4,f5,f6;
double rlx=1.0/tau;
for (n=start; n<finish; n++){
//Load data
rho_e = Den_charge[n];
rho_e = deltaT*rho_e/epsilon_LB;
rho_e = gamma*rho_e/epsilon_LB;
f0 = dist[n];
f1 = dist[2*Np+n];
@ -96,9 +98,9 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, double *Den_charge, dou
psi = f0+f2+f1+f4+f3+f6+f5;
Ex = f1-f2;
Ey = f3-f4;
Ez = f5-f6;
Ex = (f1-f2)*rlx*4.5;//NOTE the unit of electric field here is V/lu
Ey = (f3-f4)*rlx*4.5;
Ez = (f5-f6)*rlx*4.5;
ElectricField[n+0*Np] = Ex;
ElectricField[n+1*Np] = Ey;
ElectricField[n+2*Np] = Ez;

957
cpu/Stokes.cpp Normal file
View File

@ -0,0 +1,957 @@
#include <stdio.h>
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, int start, int finish, int Np)
{
double fq;
// conserved momemnts
double rho,jx,jy,jz;
double ux,uy,uz;
// non-conserved moments
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
// body force due to electric field
double rhoE;//charge density
double Ex,Ey,Ez;
// total body force
double Fx,Fy,Fz;
constexpr double mrt_V1=0.05263157894736842;
constexpr double mrt_V2=0.012531328320802;
constexpr double mrt_V3=0.04761904761904762;
constexpr double mrt_V4=0.004594820384294068;
constexpr double mrt_V5=0.01587301587301587;
constexpr double mrt_V6=0.0555555555555555555555555;
constexpr double mrt_V7=0.02777777777777778;
constexpr double mrt_V8=0.08333333333333333;
constexpr double mrt_V9=0.003341687552213868;
constexpr double mrt_V10=0.003968253968253968;
constexpr double mrt_V11=0.01388888888888889;
constexpr double mrt_V12=0.04166666666666666;
for (int n=start; n<finish; n++){
//Load data
rhoE = ChargeDensity[n];
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;
Fy = Gy + rhoE*Ey;
Fz = Gz + rhoE*Ez;
// q=0
fq = dist[n];
rho = fq;
m1 = -30.0*fq;
m2 = 12.0*fq;
// q=1
fq = dist[2*Np+n];
rho += fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jx = fq;
m4 = -4.0*fq;
m9 = 2.0*fq;
m10 = -4.0*fq;
// f2 = dist[10*Np+n];
fq = dist[1*Np+n];
rho += fq;
m1 -= 11.0*(fq);
m2 -= 4.0*(fq);
jx -= fq;
m4 += 4.0*(fq);
m9 += 2.0*(fq);
m10 -= 4.0*(fq);
// q=3
fq = dist[4*Np+n];
rho += fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jy = fq;
m6 = -4.0*fq;
m9 -= fq;
m10 += 2.0*fq;
m11 = fq;
m12 = -2.0*fq;
// q = 4
fq = dist[3*Np+n];
rho+= fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jy -= fq;
m6 += 4.0*fq;
m9 -= fq;
m10 += 2.0*fq;
m11 += fq;
m12 -= 2.0*fq;
// q=5
fq = dist[6*Np+n];
rho += fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jz = fq;
m8 = -4.0*fq;
m9 -= fq;
m10 += 2.0*fq;
m11 -= fq;
m12 += 2.0*fq;
// q = 6
fq = dist[5*Np+n];
rho+= fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jz -= fq;
m8 += 4.0*fq;
m9 -= fq;
m10 += 2.0*fq;
m11 -= fq;
m12 += 2.0*fq;
// q=7
fq = dist[8*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx += fq;
m4 += fq;
jy += fq;
m6 += fq;
m9 += fq;
m10 += fq;
m11 += fq;
m12 += fq;
m13 = fq;
m16 = fq;
m17 = -fq;
// q = 8
fq = dist[7*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx -= fq;
m4 -= fq;
jy -= fq;
m6 -= fq;
m9 += fq;
m10 += fq;
m11 += fq;
m12 += fq;
m13 += fq;
m16 -= fq;
m17 += fq;
// q=9
fq = dist[10*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx += fq;
m4 += fq;
jy -= fq;
m6 -= fq;
m9 += fq;
m10 += fq;
m11 += fq;
m12 += fq;
m13 -= fq;
m16 += fq;
m17 += fq;
// q = 10
fq = dist[9*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx -= fq;
m4 -= fq;
jy += fq;
m6 += fq;
m9 += fq;
m10 += fq;
m11 += fq;
m12 += fq;
m13 -= fq;
m16 -= fq;
m17 -= fq;
// q=11
fq = dist[12*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx += fq;
m4 += fq;
jz += fq;
m8 += fq;
m9 += fq;
m10 += fq;
m11 -= fq;
m12 -= fq;
m15 = fq;
m16 -= fq;
m18 = fq;
// q=12
fq = dist[11*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx -= fq;
m4 -= fq;
jz -= fq;
m8 -= fq;
m9 += fq;
m10 += fq;
m11 -= fq;
m12 -= fq;
m15 += fq;
m16 += fq;
m18 -= fq;
// q=13
fq = dist[14*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx += fq;
m4 += fq;
jz -= fq;
m8 -= fq;
m9 += fq;
m10 += fq;
m11 -= fq;
m12 -= fq;
m15 -= fq;
m16 -= fq;
m18 -= fq;
// q=14
fq = dist[13*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx -= fq;
m4 -= fq;
jz += fq;
m8 += fq;
m9 += fq;
m10 += fq;
m11 -= fq;
m12 -= fq;
m15 -= fq;
m16 += fq;
m18 += fq;
// q=15
fq = dist[16*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jy += fq;
m6 += fq;
jz += fq;
m8 += fq;
m9 -= 2.0*fq;
m10 -= 2.0*fq;
m14 = fq;
m17 += fq;
m18 -= fq;
// q=16
fq = dist[15*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jy -= fq;
m6 -= fq;
jz -= fq;
m8 -= fq;
m9 -= 2.0*fq;
m10 -= 2.0*fq;
m14 += fq;
m17 -= fq;
m18 += fq;
// q=17
fq = dist[18*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jy += fq;
m6 += fq;
jz -= fq;
m8 -= fq;
m9 -= 2.0*fq;
m10 -= 2.0*fq;
m14 -= fq;
m17 += fq;
m18 += fq;
// q=18
fq = dist[17*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jy -= fq;
m6 -= fq;
jz += fq;
m8 += fq;
m9 -= 2.0*fq;
m10 -= 2.0*fq;
m14 -= fq;
m17 -= fq;
m18 -= fq;
// write the velocity
ux = jx / rho;
uy = jy / rho;
uz = jz / rho;
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
//........................................................................
// READ THE DISTRIBUTIONS
// (read from opposite array due to previous swap operation)
//........................................................................
//..............incorporate external force................................................
//..............carry out relaxation process...............................................
m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho) - m2);
m4 = m4 + rlx_setB*((-0.6666666666666666*jx) - m4);
m6 = m6 + rlx_setB*((-0.6666666666666666*jy) - m6);
m8 = m8 + rlx_setB*((-0.6666666666666666*jz) - m8);
m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) - m9);
m10 = m10 + rlx_setA*(-0.5*((2*jx*jx-jy*jy-jz*jz)/rho) - m10);
m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) - m11);
m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho) - m12);
m13 = m13 + rlx_setA*((jx*jy/rho) - m13);
m14 = m14 + rlx_setA*((jy*jz/rho) - m14);
m15 = m15 + rlx_setA*((jx*jz/rho) - m15);
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
//.......................................................................................................
//.................inverse transformation......................................................
// q=0
fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2;
dist[n] = fq;
// q = 1
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10) + 0.16666666*Fx;
dist[1*Np+n] = fq;
// q=2
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10) - 0.16666666*Fx;
dist[2*Np+n] = fq;
// q = 3
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) + 0.16666666*Fy;
dist[3*Np+n] = fq;
// q = 4
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) - 0.16666666*Fy;
dist[4*Np+n] = fq;
// q = 5
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) + 0.16666666*Fz;
dist[5*Np+n] = fq;
// q = 6
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) - 0.16666666*Fz;
dist[6*Np+n] = fq;
// q = 7
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)
+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+mrt_V12*m12+0.25*m13+0.125*(m16-m17) + 0.08333333333*(Fx+Fy);
dist[7*Np+n] = fq;
// q = 8
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+mrt_V12*m12+0.25*m13+0.125*(m17-m16) - 0.08333333333*(Fx+Fy);
dist[8*Np+n] = fq;
// q = 9
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)
+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+mrt_V12*m12-0.25*m13+0.125*(m16+m17) + 0.08333333333*(Fx-Fy);
dist[9*Np+n] = fq;
// q = 10
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)
+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+mrt_V12*m12-0.25*m13-0.125*(m16+m17)- 0.08333333333*(Fx-Fy);
dist[10*Np+n] = fq;
// q = 11
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)
+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
-mrt_V12*m12+0.25*m15+0.125*(m18-m16) + 0.08333333333*(Fx+Fz);
dist[11*Np+n] = fq;
// q = 12
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)
+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
-mrt_V12*m12+0.25*m15+0.125*(m16-m18) - 0.08333333333*(Fx+Fz);
dist[12*Np+n] = fq;
// q = 13
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)
+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
-mrt_V12*m12-0.25*m15-0.125*(m16+m18) + 0.08333333333*(Fx-Fz);
dist[13*Np+n] = fq;
// q= 14
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)
+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
-mrt_V12*m12-0.25*m15+0.125*(m16+m18) - 0.08333333333*(Fx-Fz);
dist[14*Np+n] = fq;
// q = 15
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)
-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18) + 0.08333333333*(Fy+Fz);
dist[15*Np+n] = fq;
// q = 16
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)
-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17)- 0.08333333333*(Fy+Fz);
dist[16*Np+n] = fq;
// q = 17
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)
-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18) + 0.08333333333*(Fy-Fz);
dist[17*Np+n] = fq;
// q = 18
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)
-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18) - 0.08333333333*(Fy-Fz);
dist[18*Np+n] = fq;
//........................................................................
}
}
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, int start, int finish, int Np)
{
double fq;
// conserved momemnts
double rho,jx,jy,jz;
double ux,uy,uz;
// non-conserved moments
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
int nread;
// body force due to electric field
double rhoE;//charge density
double Ex,Ey,Ez;
// total body force
double Fx,Fy,Fz;
constexpr double mrt_V1=0.05263157894736842;
constexpr double mrt_V2=0.012531328320802;
constexpr double mrt_V3=0.04761904761904762;
constexpr double mrt_V4=0.004594820384294068;
constexpr double mrt_V5=0.01587301587301587;
constexpr double mrt_V6=0.0555555555555555555555555;
constexpr double mrt_V7=0.02777777777777778;
constexpr double mrt_V8=0.08333333333333333;
constexpr double mrt_V9=0.003341687552213868;
constexpr double mrt_V10=0.003968253968253968;
constexpr double mrt_V11=0.01388888888888889;
constexpr double mrt_V12=0.04166666666666666;
for (int n=start; n<finish; n++){
//Load data
rhoE = ChargeDensity[n];
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;
Fy = Gy + rhoE*Ey;
Fz = Gz + rhoE*Ez;
// q=0
fq = dist[n];
rho = fq;
m1 = -30.0*fq;
m2 = 12.0*fq;
// q=1
nread = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist)
fq = dist[nread]; // reading the f1 data into register fq
//fp = dist[10*Np+n];
rho += fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jx = fq;
m4 = -4.0*fq;
m9 = 2.0*fq;
m10 = -4.0*fq;
// f2 = dist[10*Np+n];
nread = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
fq = dist[nread]; // reading the f2 data into register fq
//fq = dist[Np+n];
rho += fq;
m1 -= 11.0*(fq);
m2 -= 4.0*(fq);
jx -= fq;
m4 += 4.0*(fq);
m9 += 2.0*(fq);
m10 -= 4.0*(fq);
// q=3
nread = neighborList[n+2*Np]; // neighbor 4
fq = dist[nread];
//fq = dist[11*Np+n];
rho += fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jy = fq;
m6 = -4.0*fq;
m9 -= fq;
m10 += 2.0*fq;
m11 = fq;
m12 = -2.0*fq;
// q = 4
nread = neighborList[n+3*Np]; // neighbor 3
fq = dist[nread];
//fq = dist[2*Np+n];
rho+= fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jy -= fq;
m6 += 4.0*fq;
m9 -= fq;
m10 += 2.0*fq;
m11 += fq;
m12 -= 2.0*fq;
// q=5
nread = neighborList[n+4*Np];
fq = dist[nread];
//fq = dist[12*Np+n];
rho += fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jz = fq;
m8 = -4.0*fq;
m9 -= fq;
m10 += 2.0*fq;
m11 -= fq;
m12 += 2.0*fq;
// q = 6
nread = neighborList[n+5*Np];
fq = dist[nread];
//fq = dist[3*Np+n];
rho+= fq;
m1 -= 11.0*fq;
m2 -= 4.0*fq;
jz -= fq;
m8 += 4.0*fq;
m9 -= fq;
m10 += 2.0*fq;
m11 -= fq;
m12 += 2.0*fq;
// q=7
nread = neighborList[n+6*Np];
fq = dist[nread];
//fq = dist[13*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx += fq;
m4 += fq;
jy += fq;
m6 += fq;
m9 += fq;
m10 += fq;
m11 += fq;
m12 += fq;
m13 = fq;
m16 = fq;
m17 = -fq;
// q = 8
nread = neighborList[n+7*Np];
fq = dist[nread];
//fq = dist[4*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx -= fq;
m4 -= fq;
jy -= fq;
m6 -= fq;
m9 += fq;
m10 += fq;
m11 += fq;
m12 += fq;
m13 += fq;
m16 -= fq;
m17 += fq;
// q=9
nread = neighborList[n+8*Np];
fq = dist[nread];
//fq = dist[14*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx += fq;
m4 += fq;
jy -= fq;
m6 -= fq;
m9 += fq;
m10 += fq;
m11 += fq;
m12 += fq;
m13 -= fq;
m16 += fq;
m17 += fq;
// q = 10
nread = neighborList[n+9*Np];
fq = dist[nread];
//fq = dist[5*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx -= fq;
m4 -= fq;
jy += fq;
m6 += fq;
m9 += fq;
m10 += fq;
m11 += fq;
m12 += fq;
m13 -= fq;
m16 -= fq;
m17 -= fq;
// q=11
nread = neighborList[n+10*Np];
fq = dist[nread];
//fq = dist[15*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx += fq;
m4 += fq;
jz += fq;
m8 += fq;
m9 += fq;
m10 += fq;
m11 -= fq;
m12 -= fq;
m15 = fq;
m16 -= fq;
m18 = fq;
// q=12
nread = neighborList[n+11*Np];
fq = dist[nread];
//fq = dist[6*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx -= fq;
m4 -= fq;
jz -= fq;
m8 -= fq;
m9 += fq;
m10 += fq;
m11 -= fq;
m12 -= fq;
m15 += fq;
m16 += fq;
m18 -= fq;
// q=13
nread = neighborList[n+12*Np];
fq = dist[nread];
//fq = dist[16*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx += fq;
m4 += fq;
jz -= fq;
m8 -= fq;
m9 += fq;
m10 += fq;
m11 -= fq;
m12 -= fq;
m15 -= fq;
m16 -= fq;
m18 -= fq;
// q=14
nread = neighborList[n+13*Np];
fq = dist[nread];
//fq = dist[7*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jx -= fq;
m4 -= fq;
jz += fq;
m8 += fq;
m9 += fq;
m10 += fq;
m11 -= fq;
m12 -= fq;
m15 -= fq;
m16 += fq;
m18 += fq;
// q=15
nread = neighborList[n+14*Np];
fq = dist[nread];
//fq = dist[17*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jy += fq;
m6 += fq;
jz += fq;
m8 += fq;
m9 -= 2.0*fq;
m10 -= 2.0*fq;
m14 = fq;
m17 += fq;
m18 -= fq;
// q=16
nread = neighborList[n+15*Np];
fq = dist[nread];
//fq = dist[8*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jy -= fq;
m6 -= fq;
jz -= fq;
m8 -= fq;
m9 -= 2.0*fq;
m10 -= 2.0*fq;
m14 += fq;
m17 -= fq;
m18 += fq;
// q=17
//fq = dist[18*Np+n];
nread = neighborList[n+16*Np];
fq = dist[nread];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jy += fq;
m6 += fq;
jz -= fq;
m8 -= fq;
m9 -= 2.0*fq;
m10 -= 2.0*fq;
m14 -= fq;
m17 += fq;
m18 += fq;
// q=18
nread = neighborList[n+17*Np];
fq = dist[nread];
//fq = dist[9*Np+n];
rho += fq;
m1 += 8.0*fq;
m2 += fq;
jy -= fq;
m6 -= fq;
jz += fq;
m8 += fq;
m9 -= 2.0*fq;
m10 -= 2.0*fq;
m14 -= fq;
m17 -= fq;
m18 -= fq;
// write the velocity
ux = jx / rho;
uy = jy / rho;
uz = jz / rho;
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
//..............incorporate external force................................................
//..............carry out relaxation process...............................................
m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho - 11*rho) - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho) - m2);
m4 = m4 + rlx_setB*((-0.6666666666666666*jx) - m4);
m6 = m6 + rlx_setB*((-0.6666666666666666*jy) - m6);
m8 = m8 + rlx_setB*((-0.6666666666666666*jz) - m8);
m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho) - m9);
m10 = m10 + rlx_setA*(-0.5*((2*jx*jx-jy*jy-jz*jz)/rho) - m10);
m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho) - m11);
m12 = m12 + rlx_setA*(-0.5*((jy*jy-jz*jz)/rho) - m12);
m13 = m13 + rlx_setA*((jx*jy/rho) - m13);
m14 = m14 + rlx_setA*((jy*jz/rho) - m14);
m15 = m15 + rlx_setA*((jx*jz/rho) - m15);
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
//.......................................................................................................
//.................inverse transformation......................................................
// q=0
fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2;
dist[n] = fq;
// q = 1
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10)+0.16666666*Fx;
nread = neighborList[n+Np];
dist[nread] = fq;
// q=2
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10) - 0.16666666*Fx;
nread = neighborList[n];
dist[nread] = fq;
// q = 3
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) + 0.16666666*Fy;
nread = neighborList[n+3*Np];
dist[nread] = fq;
// q = 4
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12) - 0.16666666*Fy;
nread = neighborList[n+2*Np];
dist[nread] = fq;
// q = 5
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) + 0.16666666*Fz;
nread = neighborList[n+5*Np];
dist[nread] = fq;
// q = 6
fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11) - 0.16666666*Fz;
nread = neighborList[n+4*Np];
dist[nread] = fq;
// q = 7
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)
+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+mrt_V12*m12+0.25*m13+0.125*(m16-m17) + 0.08333333333*(Fx+Fy);
nread = neighborList[n+7*Np];
dist[nread] = fq;
// q = 8
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+mrt_V12*m12+0.25*m13+0.125*(m17-m16) - 0.08333333333*(Fx+Fy);
nread = neighborList[n+6*Np];
dist[nread] = fq;
// q = 9
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)
+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+mrt_V12*m12-0.25*m13+0.125*(m16+m17) + 0.08333333333*(Fx-Fy);
nread = neighborList[n+9*Np];
dist[nread] = fq;
// q = 10
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)
+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+mrt_V12*m12-0.25*m13-0.125*(m16+m17)- 0.08333333333*(Fx-Fy);
nread = neighborList[n+8*Np];
dist[nread] = fq;
// q = 11
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)
+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
-mrt_V12*m12+0.25*m15+0.125*(m18-m16) + 0.08333333333*(Fx+Fz);
nread = neighborList[n+11*Np];
dist[nread] = fq;
// q = 12
fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)
+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
-mrt_V12*m12+0.25*m15+0.125*(m16-m18) - 0.08333333333*(Fx+Fz);
nread = neighborList[n+10*Np];
dist[nread]= fq;
// q = 13
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)
+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
-mrt_V12*m12-0.25*m15-0.125*(m16+m18) + 0.08333333333*(Fx-Fz);
nread = neighborList[n+13*Np];
dist[nread] = fq;
// q= 14
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)
+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
-mrt_V12*m12-0.25*m15+0.125*(m16+m18) - 0.08333333333*(Fx-Fz);
nread = neighborList[n+12*Np];
dist[nread] = fq;
// q = 15
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)
-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18) + 0.08333333333*(Fy+Fz);
nread = neighborList[n+15*Np];
dist[nread] = fq;
// q = 16
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)
-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17)- 0.08333333333*(Fy+Fz);
nread = neighborList[n+14*Np];
dist[nread] = fq;
// q = 17
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)
-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18) + 0.08333333333*(Fy-Fz);
nread = neighborList[n+17*Np];
dist[nread] = fq;
// q = 18
fq = mrt_V1*rho+mrt_V9*m1
+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)
-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18) - 0.08333333333*(Fy-Fz);
nread = neighborList[n+16*Np];
dist[nread] = fq;
}
}

190
models/IonModel.cpp
View File

@ -16,68 +16,70 @@ ScaLBL_IonModel::~ScaLBL_IonModel(){
}
void ScaLBL_IonModel::ReadParams(string filename){
//fundamental constant
kb = 1.38e-23;//Boltzmann constant;unit [J/K]
electron_charge = 1.6e-19;//electron charge;unit [C]
void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stokes,double time_conv_Stokes){
// read the input database
db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" );
ion_db = db->getDatabase( "Ions" );
// Default model parameters
//------ Load number of iteration from multiphysics controller ------//
timestepMax = num_iter;
//compute time conversion factor for ion model
time_conv = num_iter_Stokes*time_conv_Stokes/num_iter;
//-------------------------------------------------------------------//
// Universal constant
kb = 1.38e-23;//Boltzmann constant;unit [J/K]
electron_charge = 1.6e-19;//electron charge;unit [C]
//---------------------- Default model parameters --------------------------//
T = 300.0;//temperature; unit [K]
Vt = kb*T/electron_charge;//thermal voltage; unit [V]
Vt = kb*T/electron_charge;//thermal voltage; unit [Vy]
k2_inv = 4.5;//the inverse of 2nd-rank moment of D3Q7 lattice
h = 1.0;//resolution; unit: um/lu
timestepMax = 100000;
tolerance = 1.0e-8;
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]);
IonDiffusivity.push_back(1.0e-9);//user-input diffusivity has physical unit [m^2/sec]
IonValence.push_back(1);//algebraic valence charge
IonConcentration.push_back(1.0e-3);//user-input ion concentration has physical unit [mol/m^3]
//deltaT.push_back(1.0);
//tau.push_back(0.5+k2_inv*deltaT[0]*IonDiffusivisty[0]);
tau.push_back(0.5+k2_inv*time_conv/(h*1.0e-6)/(h*1.0e-6)*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( "timestepMax" )){
// timestepMax = ion_db->getScalar<int>( "timestepMax" );
//}
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" )){
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]
//NOTE: ion diffusivity has INPUT unit: [m^2/sec]
// it must be converted to LB unit: [lu^2/lt]
if (ion_db->keyExists("IonDiffusivityList")){
IonDiffusivity.clear();
IonDiffusivity = ion_db->getVector<double>( "IonDiffusivityList" );
// time relaxation parameters tau also needs update
tau.clear();
if (IonDiffusivity.size()!=number_ion_species){
ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n");
}
else{
for (int i=0; i<IonDiffusivity.size();i++){
//First, re-calculate tau
tau[i] = 0.5+k2_inv*time_conv/(h*h*1.0e-12)*IonDiffusivisty[i];
//Second, convert ion diffusivity in physical unit to LB unit
IonDiffusivity[i] = IonDiffusivity[i]*time_conv/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
}
}
}
//read ion algebric valence list
if (ion_db->keyExists("IonValenceList")){
@ -87,14 +89,34 @@ void ScaLBL_IonModel::ReadParams(string filename){
ERROR("Error: number_ion_species and IonValenceList must be the same length! \n");
}
}
//read initial ion concentration list; unit [mol/m^3]
//read initial ion concentration list; INPUT unit [mol/m^3]
//it must be converted to LB unit [mol/lu^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");
}
else{
for (int i=0; i<IonDiffusivity.size();i++){
IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
// wrong...
//if (ion_db->keyExists( "deltaT" )){
// deltaT.clear();
// tau.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");
// }
// else{//update relaxation parameter tau
// for (int i=0;i<deltaT.size();i++){
// tau[i]=0.5+k2_inv*deltaT[i]*IonDiffusivity[i];
// }
// }
//}
// Read domain parameters
if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
@ -103,15 +125,15 @@ void ScaLBL_IonModel::ReadParams(string filename){
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]);
}
if (rank==0) printf("*****************************************************\n");
if (rank==0) printf("LB Ion Transport Solver: \n");
if (rank==0) printf(" Time conversion factor: %.5g [sec/lt]\n", time_conv);
if (rank==0) printf(" Internal iteration: %i [lt]\n", timestepMax);
for (int i=0; i<number_ion_species;i++){
if (rank==0) printf(" Ion %i: LB relaxation tau = %.5g\n", i+1,tau[i]);
}
if (rank==0) printf("*****************************************************\n");
}
void ScaLBL_IonModel::SetDomain(){
@ -193,9 +215,9 @@ void ScaLBL_IonModel::ReadInput(){
}
}
// MeanFilter(Averages->SDs);
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
if (rank==0) printf("LB Ion Solver: Initialized solid phase & converting to Signed Distance function \n");
CalcDist(Distance,id_solid,*Dm);
if (rank == 0) cout << "Domain set." << endl;
if (rank == 0) cout << " Domain set." << endl;
}
void ScaLBL_IonModel::Create(){
@ -209,13 +231,13 @@ void ScaLBL_IonModel::Create(){
Mask->CommInit();
Np=Mask->PoreCount();
//...........................................................................
if (rank==0) printf ("Create ScaLBL_Communicator \n");
if (rank==0) printf ("LB Ion Solver: Create ScaLBL_Communicator \n");
// Create a communicator for the device (will use optimized layout)
// ScaLBL_Communicator ScaLBL_Comm(Mask); // original
ScaLBL_Comm = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
int Npad=(Np/16 + 2)*16;
if (rank==0) printf ("Set up memory efficient layout \n");
if (rank==0) printf ("LB Ion Solver: Set up memory efficient layout \n");
Map.resize(Nx,Ny,Nz); Map.fill(-2);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
@ -224,7 +246,7 @@ void ScaLBL_IonModel::Create(){
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
// LBM variables
if (rank==0) printf ("Allocating distributions \n");
if (rank==0) printf ("LB Ion Solver: Allocating distributions \n");
//......................device distributions.................................
int dist_mem_size = Np*sizeof(double);
int neighborSize=18*(Np*sizeof(int));
@ -235,7 +257,7 @@ void ScaLBL_IonModel::Create(){
ScaLBL_AllocateDeviceMemory((void **) &ChargeDensity, sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("Setting up device map and neighbor list \n");
if (rank==0) printf ("LB Ion Solver: Setting up device map and neighbor list \n");
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
MPI_Barrier(comm);
@ -246,10 +268,13 @@ void ScaLBL_IonModel::Initialize(){
/*
* This function initializes model
*/
if (rank==0) printf ("Initializing D3Q7 distributions for ion transport\n");
if (rank==0) printf ("LB Ion Solver: initializing D3Q7 distributions\n");
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init(&fq[ic*Np*7],&Ci[ic*Np],IonConcentration[ic],Np);
}
if (rank==0) printf ("LB Ion Solver: initializing charge density\n");
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence, number_ion_species, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence, number_ion_species, 0, ScaLBL_Comm->LastExterior(), Np);
}
void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
@ -260,13 +285,10 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
item = 1.0/item;
}
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
starttime = MPI_Wtime();
//double starttime,stoptime,cputime;
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//starttime = MPI_Wtime();
if (rank==0) printf("***************************************************\n");
if (rank==0) printf("LB-Ion Transport: timestepMax = %i\n", timestepMax);
if (rank==0) printf("***************************************************\n");
timestep=0;
while (timestep < timestepMax) {
//************************************************************************/
@ -277,22 +299,15 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_IonConcentration(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE
/* Update exterior ion concentration */
}
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++){
//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);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence, number_ion_species, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence, 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);
rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
// Set boundary conditions
@ -300,7 +315,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
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);
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
@ -309,28 +324,17 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//Update ion concentration and charge density
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FORM NORMAL
}
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
}
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);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence, number_ion_species, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence, 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);
rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
// Set boundary conditions
@ -338,25 +342,25 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
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);
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
}
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/
}
//************************************************************************/
stoptime = MPI_Wtime();
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
// Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000;
//stoptime = MPI_Wtime();
//if (rank==0) printf("-------------------------------------------------------------------\n");
//// Compute the walltime per timestep
//cputime = (stoptime - starttime)/timestep;
//// Performance obtained from each node
//double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
//if (rank==0) printf("********************************************************\n");
//if (rank==0) printf("CPU time = %f \n", cputime);
//if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
//MLUPS *= nprocs;
//if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
//if (rank==0) printf("********************************************************\n");
}

12
models/IonModel.h
View File

@ -22,24 +22,28 @@ public:
~ScaLBL_IonModel();
// functions in they should be run
void ReadParams(string filename);
void ReadParams(string filename,int num_iter,int num_iter_Stokes,double time_conv_Stokes);
void ReadParams(std::shared_ptr<Database> db0);
void SetDomain();
void ReadInput();
void Create();
void Initialize();
void Run(double *Velocity);
void VelocityField();
void Run(double *Velocity, double *ElectricField);
bool Restart,pBC;
int timestep,timestepMax;
int BoundaryCondition;
double h;//domain resolution, unit [um/lu]
double time_conv;
double kb,electron_charge,T,Vt;
double k2_inv;
double tolerance;
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> deltaT;
vector<double> tau;
int Nx,Ny,Nz,N,Np;

57
models/MultiPhysController.cpp Normal file
View File

@ -0,0 +1,57 @@
#include "models/MultiPhysController.h"
ScaLBL_Multiphys_Controller::ScaLBL_Multiphys_Controller(int RANK, int NP, MPI_Comm COMM):
rank(RANK),Restart(0),timestepMax(0),num_iter_Stokes(0),num_iter_Ion(0),SchmidtNum(0),comm(COMM)
{
}
ScaLBL_Multiphys_Controller::~ScaLBL_Multiphys_Controller(){
}
void ScaLBL_Multiphys_Controller::ReadParams(string filename){
// read the input database
db = std::make_shared<Database>( filename );
study_db = db->getDatabase( "MultiphysController" );
// Default parameters
timestepMax = 10000;
Restart = false;
SchmidtNum = 1.0;
num_iter_Stokes=1;
num_iter_Ion=1;
// load input parameters
if (study_db->keyExists( "timestepMax" )){
timestepMax = study_db->getScalar<int>( "timestepMax" );
}
if (study_db->keyExists( "Schmidt_Number" )){
SchmidtNum = study_db->getScalar<double>( "Schmidt_Number" );
}
// recalculate relevant parameters
if (SchmidtNum>1){
num_iter_Stokes = int(round(SchmidtNum/2)*2);
num_iter_Ion = 1;
}
else if (SchmidtNum>0 && SchmidtNum<1){
num_iter_Ion = int(round((1.0/SchmidtNum)/2)*2);
num_iter_Stokes = 1;
}
else{
ERROR("Error: SchmidtNum (Schmidt number) must be a positive number! \n");
}
// load input parameters
// in case user wants to have an absolute control over the iternal iteration
if (study_db->keyExists( "num_iter_Ion" )){
num_iter_Ion = study_db->getScalar<int>( "num_iter_Ion" );
}
if (study_db->keyExists( "num_iter_Stokes" )){
num_iter_Stokes = study_db->getScalar<int>( "num_iter_Stokes" );
}
}

47
models/MultiPhysController.h Normal file
View File

@ -0,0 +1,47 @@
/*
* Multiphysics controller that coordinates the coupling between different models
*/
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
class ScaLBL_Multiphys_Controller{
public:
ScaLBL_Multiphys_Controller(int RANK, int NP, MPI_Comm COMM);
~ScaLBL_Multiphys_Controller();
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
bool Restart;
//int timestep;
int timestepMax;
int num_iter_Stokes;
int num_iter_Ion;
double SchmidtNum;//Schmidt number = kinematic_viscosity/mass_diffusivity
// input database
std::shared_ptr<Database> db;
std::shared_ptr<Database> study_db;
private:
MPI_Comm comm;
// filenames
char LocalRankString[8];
char LocalRankFilename[40];
char LocalRestartFile[40];
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
};

81
models/PoissonSolver.cpp
View File

@ -23,8 +23,8 @@ void ScaLBL_Poisson::ReadParams(string filename){
electric_db = db->getDatabase( "Poisson" );
k2_inv = 4.5;//the inverse of 2nd-rank moment of D3Q7 lattice
deltaT = 0.3;//time step of LB-Poisson equation
tau = 0.5+k2_inv*deltaT;
gamma = 0.3;//time step of LB-Poisson equation
tau = 0.5+k2_inv*gamma;
timestepMax = 100000;
tolerance = 1.0e-6;//stopping criterion for obtaining steady-state electricla potential
h = 1.0;//resolution; unit: um/lu
@ -44,8 +44,8 @@ void ScaLBL_Poisson::ReadParams(string filename){
if (electric_db->keyExists( "tolerance" )){
tolerance = electric_db->getScalar<double>( "tolerance" );
}
if (electric_db->keyExists( "deltaT" )){
deltaT = electric_db->getScalar<double>( "deltaT" );
if (electric_db->keyExists( "gamma" )){
gamma = electric_db->getScalar<double>( "gamma" );
}
if (electric_db->keyExists( "epsilonR" )){
epsilonR = electric_db->getScalar<double>( "epsilonR" );
@ -60,8 +60,12 @@ void ScaLBL_Poisson::ReadParams(string filename){
//Re-calcualte model parameters if user updates input
epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)]
epsilon_LB = epsilon0_LB*epsilonR;//electrical permittivity
tau = 0.5+k2_inv*deltaT;
tau = 0.5+k2_inv*gamma;
if (rank==0) printf("***********************************************************************************\n");
if (rank==0) printf("LB-Poisson Solver: steady-state MaxTimeStep = %i; steady-state tolerance = %.3g \n", timestepMax,tolerance);
if (rank==0) printf(" LB relaxation tau = %.5g \n", tau);
if (rank==0) printf("***********************************************************************************\n");
}
void ScaLBL_Poisson::SetDomain(){
Dm = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // full domain for analysis
@ -143,9 +147,9 @@ void ScaLBL_Poisson::ReadInput(){
}
}
// MeanFilter(Averages->SDs);
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
if (rank==0) printf("LB-Poisson Solver: Initialized solid phase & converting to Signed Distance function \n");
CalcDist(Distance,id_solid,*Dm);
if (rank == 0) cout << "Domain set." << endl;
if (rank == 0) cout << " Domain set." << endl;
}
void ScaLBL_Poisson::Create(){
@ -159,13 +163,13 @@ void ScaLBL_Poisson::Create(){
Mask->CommInit();
Np=Mask->PoreCount();
//...........................................................................
if (rank==0) printf ("Create ScaLBL_Communicator \n");
if (rank==0) printf ("LB-Poisson Solver: Create ScaLBL_Communicator \n");
// Create a communicator for the device (will use optimized layout)
// ScaLBL_Communicator ScaLBL_Comm(Mask); // original
ScaLBL_Comm = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
int Npad=(Np/16 + 2)*16;
if (rank==0) printf ("Set up memory efficient layout \n");
if (rank==0) printf ("LB-Poisson Solver: Set up memory efficient layout \n");
Map.resize(Nx,Ny,Nz); Map.fill(-2);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
@ -174,7 +178,7 @@ void ScaLBL_Poisson::Create(){
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
// LBM variables
if (rank==0) printf ("Allocating distributions \n");
if (rank==0) printf ("LB-Poisson Solver: Allocating distributions \n");
//......................device distributions.................................
int dist_mem_size = Np*sizeof(double);
int neighborSize=18*(Np*sizeof(int));
@ -185,7 +189,7 @@ void ScaLBL_Poisson::Create(){
ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("Setting up device map and neighbor list \n");
if (rank==0) printf ("LB-Poisson Solver: Setting up device map and neighbor list \n");
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
MPI_Barrier(comm);
@ -196,22 +200,17 @@ void ScaLBL_Poisson::Initialize(){
/*
* This function initializes model
*/
if (rank==0) printf ("Initializing D3Q7 distributions for LB-Poisson solver\n");
if (rank==0) printf ("LB-Poisson Solver: initializing D3Q7 distributions\n");
ScaLBL_D3Q7_Poisson_Init(fq, Np);
}
void ScaLBL_Poisson::Run(double *ChargeDensity){
//LB-related parameter
double rlx = 1.0/tau;
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
starttime = MPI_Wtime();
if (rank==0) printf("***************************************************************************\n");
if (rank==0) printf("LB-Poisson Solver: timestepMax = %i; steady-state tolerance = %.3g \n", timestepMax,tolerance);
if (rank==0) printf("***************************************************************************\n");
//double starttime,stoptime,cputime;
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//starttime = MPI_Wtime();
timestep=0;
double error = 1.0;
double psi_avg_previous = 0.0;
@ -220,21 +219,21 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
// *************ODD TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, rlx, epsilon_LB, deltaT, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(fq, 0); //WRITE INTO OPPOSITE
// Set boundary conditions
/* ... */
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, rlx, epsilon_LB, deltaT, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// *************EVEN TIMESTEP*************//
timestep++;
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FORM NORMAL
ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, rlx, epsilon_LB, deltaT, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(fq, 0); //WRITE INTO OPPOSITE
// Set boundary conditions
/* ... */
ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, rlx, epsilon_LB, deltaT, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, gamma, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/
@ -267,24 +266,22 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
psi_avg_previous = psi_avg;
}
}
//************************************************************************/
stoptime = MPI_Wtime();
if (rank==0) printf("LB-Poission Solver: a steady-state solution is obtained\n");
if (rank==0) printf("---------------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
// Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("******************* LB-Poisson Solver Statistics ********************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("*********************************************************************\n");
//************************************************************************/
//stoptime = MPI_Wtime();
////if (rank==0) printf("LB-Poission Solver: a steady-state solution is obtained\n");
////if (rank==0) printf("---------------------------------------------------------------------------\n");
//// Compute the walltime per timestep
//cputime = (stoptime - starttime)/timestep;
//// Performance obtained from each node
//double MLUPS = double(Np)/cputime/1000000;
//if (rank==0) printf("******************* LB-Poisson Solver Statistics ********************\n");
//if (rank==0) printf("CPU time = %f \n", cputime);
//if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
//MLUPS *= nprocs;
//if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
//if (rank==0) printf("*********************************************************************\n");
}
//void ScaLBL_Poisson::get_ElectricField(){
//// ???
//}

543
models/StokesModel.cpp Normal file
View File

@ -0,0 +1,543 @@
/*
* Multi-relaxation time LBM Model
*/
#include "models/StokesModel.h"
#include "analysis/distance.h"
#include "common/ReadMicroCT.h"
ScaLBL_StokesModel::ScaLBL_StokesModel(int RANK, int NP, MPI_Comm COMM):
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),mu(0),
Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{
}
ScaLBL_StokesModel::~ScaLBL_StokesModel(){
}
void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
// read the input database
db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" );
stokes_db = db->getDatabase( "Stokes" );
//------ Load number of iteration from multiphysics controller ------//
timestepMax = num_iter;
//-------------------------------------------------------------------//
//---------------------- Default model parameters --------------------------//
nu_phys = 1.004e-6;//by default use water kinematic viscosity at 20C; unit [m^2/sec]
h = 1.0;//image resolution;[um]
tau = 1.0;
mu = (tau-0.5)/3.0;//LB kinematic viscosity;unit [lu^2/lt]
time_conv = h*h*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
tolerance = 1.0e-8;
Fx = Fy = 0.0;
Fz = 1.0e-5;
//--------------------------------------------------------------------------//
// Single-fluid Navier-Stokes Model parameters
//if (stokes_db->keyExists( "timestepMax" )){
// timestepMax = stokes_db->getScalar<int>( "timestepMax" );
//}
if (stokes_db->keyExists( "tolerance" )){
tolerance = stokes_db->getScalar<double>( "tolerance" );
}
if (stokes_db->keyExists( "tau" )){
tau = stokes_db->getScalar<double>( "tau" );
}
if (stokes_db->keyExists( "nu_phys" )){
nu_phys = stokes_db->getScalar<double>( "nu_phys" );
}
if (stokes_db->keyExists( "F" )){
Fx = stokes_db->getVector<double>( "F" )[0];
Fy = stokes_db->getVector<double>( "F" )[1];
Fz = stokes_db->getVector<double>( "F" )[2];
}
if (stokes_db->keyExists( "Restart" )){
Restart = stokes_db->getScalar<bool>( "Restart" );
}
if (stokes_db->keyExists( "din" )){
din = stokes_db->getScalar<double>( "din" );
}
if (stokes_db->keyExists( "dout" )){
dout = stokes_db->getScalar<double>( "dout" );
}
if (stokes_db->keyExists( "flux" )){
flux = stokes_db->getScalar<double>( "flux" );
}
// Read domain parameters
if (domain_db->keyExists( "BC" )){
BoundaryCondition = domain_db->getScalar<int>( "BC" );
}
if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
h = domain_db->getScalar<double>( "voxel_length" );
}
// Re-calculate model parameters due to parameter read
mu=(tau-0.5)/3.0;
time_conv = h*h*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
if (rank==0) printf("*****************************************************\n");
if (rank==0) printf("LB Single-Fluid Navier-Stokes Solver: \n");
if (rank==0) printf(" Time conversion factor: %.5g [sec/lt]\n", time_conv);
if (rank==0) printf(" Internal iteration: %i [lt]\n", timestepMax);
if (rank==0) printf("*****************************************************\n");
}
void ScaLBL_StokesModel::SetDomain(){
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
// domain parameters
Nx = Dm->Nx;
Ny = Dm->Ny;
Nz = Dm->Nz;
Lx = Dm->Lx;
Ly = Dm->Ly;
Lz = Dm->Lz;
N = Nx*Ny*Nz;
Distance.resize(Nx,Ny,Nz);
Velocity_x.resize(Nx,Ny,Nz);
Velocity_y.resize(Nx,Ny,Nz);
Velocity_z.resize(Nx,Ny,Nz);
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
MPI_Barrier(comm);
Dm->CommInit();
MPI_Barrier(comm);
rank = Dm->rank();
nprocx = Dm->nprocx();
nprocy = Dm->nprocy();
nprocz = Dm->nprocz();
}
void ScaLBL_StokesModel::ReadInput(){
sprintf(LocalRankString,"%05d",Dm->rank());
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
if (domain_db->keyExists( "Filename" )){
auto Filename = domain_db->getScalar<std::string>( "Filename" );
Mask->Decomp(Filename);
}
else if (domain_db->keyExists( "GridFile" )){
// Read the local domain data
auto input_id = readMicroCT( *domain_db, comm );
// Fill the halo (assuming GCW of 1)
array<int,3> size0 = { (int) input_id.size(0), (int) input_id.size(1), (int) input_id.size(2) };
ArraySize size1 = { (size_t) Mask->Nx, (size_t) Mask->Ny, (size_t) Mask->Nz };
ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
fillHalo<signed char> fill( comm, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
Array<signed char> id_view;
id_view.viewRaw( size1, Mask->id );
fill.copy( input_id, id_view );
fill.fill( id_view );
}
else{
Mask->ReadIDs();
}
// Generate the signed distance map
// Initialize the domain and communication
Array<char> id_solid(Nx,Ny,Nz);
// Solve for the position of the solid phase
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
// Initialize the solid phase
if (Mask->id[n] > 0) id_solid(i,j,k) = 1;
else id_solid(i,j,k) = 0;
}
}
}
// Initialize the signed distance function
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
// Initialize distance to +/- 1
Distance(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
}
}
}
// MeanFilter(Averages->SDs);
if (rank==0) printf("LB Single-Fluid Solver: initialized solid phase & converting to Signed Distance function \n");
CalcDist(Distance,id_solid,*Dm);
if (rank == 0) cout << " Domain set." << endl;
}
void ScaLBL_StokesModel::Create(){
/*
* This function creates the variables needed to run a LBM
*/
int rank=Mask->rank();
//.........................................................
// Initialize communication structures in averaging domain
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
Mask->CommInit();
Np=Mask->PoreCount();
//...........................................................................
if (rank==0) printf ("LB Single-Fluid Solver: Create ScaLBL_Communicator \n");
// Create a communicator for the device (will use optimized layout)
// ScaLBL_Communicator ScaLBL_Comm(Mask); // original
ScaLBL_Comm = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
int Npad=(Np/16 + 2)*16;
if (rank==0) printf ("LB Single-Fluid Solver: Set up memory efficient layout \n");
Map.resize(Nx,Ny,Nz); Map.fill(-2);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
MPI_Barrier(comm);
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
// LBM variables
if (rank==0) printf ("LB Single-Fluid Solver: Allocating distributions \n");
//......................device distributions.................................
int dist_mem_size = Np*sizeof(double);
int neighborSize=18*(Np*sizeof(int));
//...........................................................................
ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("LB Single-Fluid Solver: Setting up device map and neighbor list \n");
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
MPI_Barrier(comm);
}
void ScaLBL_StokesModel::Initialize(){
/*
* This function initializes model
*/
if (rank==0) printf("LB Single-Fluid Solver: Initializing distributions \n");
if (rank==0) printf("****************************************************************\n");
ScaLBL_D3Q19_Init(fq, Np);
}
void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
double rlx_setA=1.0/tau;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
timestep = 0;
while (timestep < timestepMax) {
//************************************************************************/
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
if (BoundaryCondition == 3){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 4){
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 5){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
if (BoundaryCondition == 3){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 4){
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 5){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/
}
}
//void ScaLBL_StokesModel::computeVelocity_phys(){
// ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
// ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//}
void ScaLBL_StokesModel::Run(){
double rlx_setA=1.0/tau;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
Minkowski Morphology(Mask);
if (rank==0){
bool WriteHeader=false;
FILE *log_file = fopen("Permeability.csv","r");
if (log_file != NULL)
fclose(log_file);
else
WriteHeader=true;
if (WriteHeader){
log_file = fopen("Permeability.csv","a+");
fprintf(log_file,"time Fx Fy Fz mu Vs As Js Xs vx vy vz k\n");
fclose(log_file);
}
}
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
starttime = MPI_Wtime();
if (rank==0) printf("****************************************************************\n");
if (rank==0) printf("LB Single-Fluid Navier-Stokes Solver: timestepMax = %i\n", timestepMax);
if (rank==0) printf("****************************************************************\n");
timestep=0;
double error = 1.0;
double flow_rate_previous = 0.0;
while (timestep < timestepMax && error > tolerance) {
//************************************************************************/
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
if (BoundaryCondition == 3){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 4){
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 5){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_MRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
if (BoundaryCondition == 3){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 4){
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 5){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAeven_MRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/
if (timestep%1000==0){
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Velocity_x);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Velocity_y);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Velocity_z);
double count_loc=0;
double count;
double vax,vay,vaz;
double vax_loc,vay_loc,vaz_loc;
vax_loc = vay_loc = vaz_loc = 0.f;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
if (Distance(i,j,k) > 0){
vax_loc += Velocity_x(i,j,k);
vay_loc += Velocity_y(i,j,k);
vaz_loc += Velocity_z(i,j,k);
count_loc+=1.0;
}
}
}
}
MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
vax /= count;
vay /= count;
vaz /= count;
double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
double dir_x = Fx/force_mag;
double dir_y = Fy/force_mag;
double dir_z = Fz/force_mag;
if (force_mag == 0.0){
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
dir_z = 1.0;
force_mag = 1.0;
}
double flow_rate = (vax*dir_x + vay*dir_y + vaz*dir_z);
error = fabs(flow_rate - flow_rate_previous) / fabs(flow_rate);
flow_rate_previous = flow_rate;
//if (rank==0) printf("Computing Minkowski functionals \n");
Morphology.ComputeScalar(Distance,0.f);
//Morphology.PrintAll();
double mu = (tau-0.5)/3.f;
double Vs = Morphology.V();
double As = Morphology.A();
double Hs = Morphology.H();
double Xs = Morphology.X();
Vs=sumReduce( Dm->Comm, Vs);
As=sumReduce( Dm->Comm, As);
Hs=sumReduce( Dm->Comm, Hs);
Xs=sumReduce( Dm->Comm, Xs);
double h = Dm->voxel_length;
double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag;
if (rank==0) {
printf(" %f\n",absperm);
FILE * log_file = fopen("Permeability.csv","a");
fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu,
h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm);
fclose(log_file);
}
}
}
//************************************************************************/
stoptime = MPI_Wtime();
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
// Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
}
void ScaLBL_StokesModel::VelocityField(){
/* Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],3*SIZE);
memcpy(Morphology.SDn.data(), Distance.data(), Nx*Ny*Nz*sizeof(double));
Morphology.Initialize();
Morphology.UpdateMeshValues();
Morphology.ComputeLocal();
Morphology.Reduce();
double count_loc=0;
double count;
double vax,vay,vaz;
double vax_loc,vay_loc,vaz_loc;
vax_loc = vay_loc = vaz_loc = 0.f;
for (int n=0; n<ScaLBL_Comm->LastExterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
count_loc+=1.0;
}
for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
count_loc+=1.0;
}
MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
vax /= count;
vay /= count;
vaz /= count;
double mu = (tau-0.5)/3.f;
if (rank==0) printf("Fx Fy Fz mu Vs As Js Xs vx vy vz\n");
if (rank==0) printf("%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",Fx, Fy, Fz, mu,
Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz);
*/
std::vector<IO::MeshDataStruct> visData;
fillHalo<double> fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1);
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
auto SignDistVar = std::make_shared<IO::Variable>();
IO::initialize("","silo","false");
// Create the MeshDataStruct
visData.resize(1);
visData[0].meshName = "domain";
visData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VariableType::VolumeVariable;
SignDistVar->dim = 1;
SignDistVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(SignDistVar);
VxVar->name = "Velocity_x";
VxVar->type = IO::VariableType::VolumeVariable;
VxVar->dim = 1;
VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VxVar);
VyVar->name = "Velocity_y";
VyVar->type = IO::VariableType::VolumeVariable;
VyVar->dim = 1;
VyVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VyVar);
VzVar->name = "Velocity_z";
VzVar->type = IO::VariableType::VolumeVariable;
VzVar->dim = 1;
VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VzVar);
Array<double>& SignData = visData[0].vars[0]->data;
Array<double>& VelxData = visData[0].vars[1]->data;
Array<double>& VelyData = visData[0].vars[2]->data;
Array<double>& VelzData = visData[0].vars[3]->data;
ASSERT(visData[0].vars[0]->name=="SignDist");
ASSERT(visData[0].vars[1]->name=="Velocity_x");
ASSERT(visData[0].vars[2]->name=="Velocity_y");
ASSERT(visData[0].vars[3]->name=="Velocity_z");
fillData.copy(Distance,SignData);
fillData.copy(Velocity_x,VelxData);
fillData.copy(Velocity_y,VelyData);
fillData.copy(Velocity_z,VelzData);
IO::writeData( timestep, visData, Dm->Comm );
}

75
models/StokesModel.h Normal file
View File

@ -0,0 +1,75 @@
/*
* Multi-relaxation time LBM Model
*/
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
class ScaLBL_StokesModel{
public:
ScaLBL_StokesModel(int RANK, int NP, MPI_Comm COMM);
~ScaLBL_StokesModel();
// functions in they should be run
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
void SetDomain();
void ReadInput();
void Create();
void Initialize();
void Run();
void VelocityField();
bool Restart,pBC;
int timestep,timestepMax;
int BoundaryCondition;
double tau,mu;
double Fx,Fy,Fz,flux;
double din,dout;
double tolerance;
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
double Lx,Ly,Lz;
std::shared_ptr<Domain> Dm; // this domain is for analysis
std::shared_ptr<Domain> Mask; // this domain is for lbm
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
// input database
std::shared_ptr<Database> db;
std::shared_ptr<Database> domain_db;
std::shared_ptr<Database> stokes_db;
IntArray Map;
DoubleArray Distance;
int *NeighborList;
double *fq;
double *Velocity;
double *Pressure;
//Minkowski Morphology;
DoubleArray Velocity_x;
DoubleArray Velocity_y;
DoubleArray Velocity_z;
private:
MPI_Comm comm;
// filenames
char LocalRankString[8];
char LocalRankFilename[40];
char LocalRestartFile[40];
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
};

137
tests/lbpm_electrokinetic_dfh_simulator.cpp
View File

@ -5,74 +5,99 @@
#include <exception>
#include <stdexcept>
#include <fstream>
#include <math.h>
#include "models/DFHModel.h"
#include "models/StokesModel.h"
#include "models/IonModel.h"
#include "models/PoissonSolver.h"
//#define WRE_SURFACES
/*
* Simulator for two-phase flow in porous media
* James E. McClure 2013-2014
*/
#include "models/MultiPhysController.h"
using namespace std;
//*************************************************************************
// Implementation of Two-Phase Immiscible LBM using CUDA
//*************************************************************************
//***************************************************************************
// Implementation of Multiphysics simulator using lattice-Boltzmann method
//***************************************************************************
int main(int argc, char **argv)
{
// Initialize MPI
int provided_thread_support = -1;
MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
MPI_Comm comm;
MPI_Comm_dup(MPI_COMM_WORLD,&comm);
int rank = comm_rank(comm);
int nprocs = comm_size(comm);
if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE )
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
// Initialize MPI
int provided_thread_support = -1;
MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
MPI_Comm comm;
MPI_Comm_dup(MPI_COMM_WORLD,&comm);
int rank = comm_rank(comm);
int nprocs = comm_size(comm);
if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
}
if (rank == 0){
printf("********************************************************\n");
printf("Running Electrokinetic LBM Simulator \n");
printf("********************************************************\n");
}
PROFILE_ENABLE(1);
//PROFILE_ENABLE_TRACE();
//PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START("Main");
Utilities::setErrorHandlers();
// Limit scope so variables that contain communicators will free before MPI_Finialize
{
if (rank == 0){
printf("********************************************************\n");
printf("Running Electrokinetic LBM Simulator \n");
printf("********************************************************\n");
}
//PROFILE_ENABLE_TRACE();
//PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START("Main");
Utilities::setErrorHandlers();
auto filename = argv[1];
ScaLBL_DFHModel DFHModel(rank,nprocs,comm);
ScaLBL_IonModel IonModel(rank,nprocs,comm);
ScaLBL_Poisson PoissonSolver(rank,nprocs,comm);
DFHModel.ReadParams(filename);
DFHModel.SetDomain();
DFHModel.ReadInput();
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.Run(IonModel.ChargeDensity); // Solve the N-S equations to get velocity
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);
auto filename = argv[1];
ScaLBL_StokesModel StokesModel(rank,nprocs,comm);
ScaLBL_IonModel IonModel(rank,nprocs,comm);
ScaLBL_Poisson PoissonSolver(rank,nprocs,comm);
ScaLBL_Multiphys_Controller Study(rank,nprocs,comm);//multiphysics controller coordinating multi-model coupling
DFHModel.WriteDebug();
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_color_simulator",1);
// ****************************************************
MPI_Barrier(comm);
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Comm_free(&comm);
MPI_Finalize();
// Load controller information
Study.ReadParams(filename);
// Initialize LB Navier-Stokes model
StokesModel.ReadParams(filename,Study.num_iter_Stokes);
StokesModel.SetDomain();
StokesModel.ReadInput();
StokesModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables
StokesModel.Initialize(); // initializing the model will set initial conditions for variables
// Initialize LB-Ion model
IonModel.ReadParams(filename,Study.num_iter_Ion,Study.num_iter_Stokes,StokesModel.time_conv);
IonModel.SetDomain();
IonModel.ReadInput();
IonModel.Create();
IonModel.Initialize();
// Initialize LB-Poisson model
PoissonSolver.ReadParams(filename);
PoissonSolver.SetDomain();
PoissonSolver.ReadInput();
PoissonSolver.Create();
PoissonSolver.Initialize();
int timestep=0;
while (timestep < Study.timestepMax){
timestep++;
PoissonSolver.Run(IonModel.ChargeDensity);//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
//--------------------------------------------
//potentially leave analysis module for future
//--------------------------------------------
}
//StokesModel.WriteDebug();
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_electrokinetic_simulator",1);
// ****************************************************
MPI_Barrier(comm);
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Comm_free(&comm);
MPI_Finalize();
}