Merge branch 'master' into SYCL

This commit is contained in:
James McClure 2023-04-03 15:49:34 -04:00
commit e7b9407586
30 changed files with 1497 additions and 132 deletions

View File

@ -430,11 +430,11 @@ void SubPhase::Basic() {
//double total_flow_rate = water_flow_rate + not_water_flow_rate;
//double fractional_flow = water_flow_rate / total_flow_rate;
double h = Dm->voxel_length;
double krn = h * h * nu_n * Porosity * not_water_flow_rate / force_mag;
double krw = h * h * nu_w * Porosity* water_flow_rate / force_mag;
double krn = h * h * nu_n * Porosity* Porosity * not_water_flow_rate / force_mag;
double krw = h * h * nu_w * Porosity* Porosity* water_flow_rate / force_mag;
/* not counting films */
double krnf = krn - h * h * nu_n * Porosity * not_water_film_flow_rate / force_mag;
double krwf = krw - h * h * nu_w * Porosity * water_film_flow_rate / force_mag;
double krnf = krn - h * h * nu_n * Porosity* Porosity * not_water_film_flow_rate / force_mag;
double krwf = krw - h * h * nu_w * Porosity* Porosity * water_film_flow_rate / force_mag;
double eff_pressure = 1.0 / (krn + krw); // effective pressure drop
fprintf(TIMELOG,

View File

@ -217,18 +217,25 @@ void Domain::read_swc(const std::string &Filename) {
double z = k*voxel_length + start_z;
double distance;
double s = ((x-xp)*alpha+(y-yp)*beta+(z-zp)*gamma) / (alpha*alpha + beta*beta + gamma*gamma);
if (s > length){
distance = ri - sqrt((x-xi)*(x-xi) + (y-yi)*(y-yi) + (z-zi)*(z-zi));
double s = ((x-xp)*alpha+(y-yp)*beta+(z-zp)*gamma) / (alpha*alpha + beta*beta + gamma*gamma);
double di = ri - sqrt((x-xi)*(x-xi) + (y-yi)*(y-yi) + (z-zi)*(z-zi));
double dp = rp - sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp) + (z-zp)*(z-zp));
if (s > length ){
distance = di;
}
else if (s < 0.0){
distance = rp - sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp) + (z-zp)*(z-zp));
distance = dp;
}
else {
// linear variation for radius
double radius = rp + (ri - rp)*s/length;
distance = radius - sqrt((x-xp-alpha*s)*(x-xp-alpha*s) + (y-yp-beta*s)*(y-yp-beta*s) + (z-zp-gamma*s)*(z-zp-gamma*s));
}
if (distance < di) distance = di;
if (distance < dp) distance = dp;
if ( distance > 0.0 ){
/* label the voxel */
//id[k*Nx*Ny + j*Nx + i] = label;

View File

@ -667,6 +667,74 @@ int Membrane::Create(DoubleArray &Distance, IntArray &Map){
return mlink;
}
void Membrane::Write(string filename){
int mlink = membraneLinkCount;
std::ofstream ofs (filename, std::ofstream::out);
/* Create local copies of membrane data structures */
double *tmpMembraneCoef; // mass transport coefficient for the membrane
tmpMembraneCoef = new double [2*mlink*sizeof(double)];
ScaLBL_CopyToHost(tmpMembraneCoef, MembraneCoef, 2*mlink*sizeof(double));
int i,j,k;
for (int m=0; m<mlink; m++){
double a1 = tmpMembraneCoef[2*m];
double a2 = tmpMembraneCoef[2*m+1];
int m1 = membraneLinks[2*m]%Np;
int m2 = membraneLinks[2*m+1]%Np;
// map index to global i,j,k
k = m1/(Nx*Ny); j = (m1-Nx*Ny*k)/Nx; i = m1-Nx*Ny*k-Nx*j;
ofs << i << " " << j << " " << k << " "<< a1 ;
k = m2/(Nx*Ny); j = (m2-Nx*Ny*k)/Nx; i = m2-Nx*Ny*k-Nx*j;
ofs << i << " " << j << " " << k << " "<< a2 << endl;
}
ofs.close();
/*FILE *VELX_FILE;
sprintf(LocalRankFilename, "Velocity_X.%05i.raw", rank);
VELX_FILE = fopen(LocalRankFilename, "wb");
fwrite(PhaseField.data(), 8, N, VELX_FILE);
fclose(VELX_FILE);
*/
delete [] tmpMembraneCoef;
}
void Membrane::Read(string filename){
int mlink = membraneLinkCount;
/* Create local copies of membrane data structures */
double *tmpMembraneCoef; // mass transport coefficient for the membrane
tmpMembraneCoef = new double [2*mlink*sizeof(double)];
FILE *fid = fopen(filename.c_str(), "r");
INSIST(fid != NULL, "Error opening membrane file \n");
//........read the spheres..................
// We will read until a blank like or end-of-file is reached
int count = 0;
int i,j,k;
int ii,jj,kk;
double a1, a2;
while (fscanf(fid, "%i,%i,%i,%lf,%i,%i,%i,%lf,\n", &i, &j, &k, &a1, &ii, &jj, &kk, &a2) == 8){
printf("%i, %i, %i, %lf \n", i,j,k, a2);
count++;
}
if (count != mlink){
printf("WARNING (Membrane::Read): number of file lines does not match number of links \n");
}
fclose(fid);
ScaLBL_CopyToDevice(MembraneCoef, tmpMembraneCoef, 2*mlink*sizeof(double));
/*FILE *VELX_FILE;
sprintf(LocalRankFilename, "Velocity_X.%05i.raw", rank);
VELX_FILE = fopen(LocalRankFilename, "wb");
fwrite(PhaseField.data(), 8, N, VELX_FILE);
fclose(VELX_FILE);
*/
delete [] tmpMembraneCoef;
}
int Membrane::D3Q7_MapRecv(int Cqx, int Cqy, int Cqz, int *d3q19_recvlist,
int count, int *membraneRecvLabels, DoubleArray &Distance, int *dvcMap){

View File

@ -92,6 +92,18 @@ public:
* @param Map - mapping between regular layout and compact layout
*/
int Create(DoubleArray &Distance, IntArray &Map);
/**
* \brief Write membrane data to output file
* @param filename - name of file to save
*/
void Write(string filename);
/**
* \brief Read membrane data from input file
* @param filename - name of file to save
*/
void Read(string filename);
void SendD3Q7AA(double *dist);
void RecvD3Q7AA(double *dist);

View File

@ -2472,6 +2472,29 @@ void ScaLBL_Communicator::D3Q7_Poisson_Potential_BC_Z(int *neighborList, double
}
}
void ScaLBL_Communicator::D3Q19_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time){
if (kproc == 0) {
if (time%2==0){
ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(dvcSendList_z, fq, Vin, sendCount_z, N);
}
else{
ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(neighborList, dvcSendList_z, fq, Vin, sendCount_z, N);
}
}
}
void ScaLBL_Communicator::D3Q19_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time){
if (kproc == nprocz-1){
if (time%2==0){
ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(dvcSendList_Z, fq, Vout, sendCount_Z, N);
}
else{
ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(neighborList, dvcSendList_Z, fq, Vout, sendCount_Z, N);
}
}
}
void ScaLBL_Communicator::Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin){
if (kproc == 0) {
ScaLBL_Poisson_D3Q7_BC_z(dvcSendList_z, Map, Psi, Vin, sendCount_z);

View File

@ -396,6 +396,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist,
extern "C" void ScaLBL_D3Q19_Poisson_getElectricField(double *dist, double *ElectricField, double tau, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vout, int count, int Np);
// LBM Stokes Model (adapted from MRT model)
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np);
@ -806,6 +812,8 @@ public:
double D3Q19_Flux_BC_z(int *neighborList, double *fq, double flux, int time);
void D3Q7_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time);
void D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time);
void D3Q19_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time);
void D3Q19_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time);
void Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin);
void Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout);
void D3Q7_Ion_Concentration_BC_z(int *neighborList, double *fq, double Cin, int time);

View File

@ -181,6 +181,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_z(int *d_neighborList,
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList,
int *list,
double *dist,
@ -213,6 +214,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList,
}
}
extern "C" void ScaLBL_Poisson_D3Q7_BC_z(int *list, int *Map, double *Psi,
double Vin, int count) {
int idx, n, nm;

View File

@ -631,6 +631,8 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson_ElectricPotential(
}
}
extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map,
double *dist, double *Den_charge,
double *Psi, double *ElectricField,
@ -915,6 +917,92 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist,
}
}
extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int nread, nr5;
for (int idx = 0; idx < count; idx++) {
int n = list[idx];
dist[6 * Np + n] = W1*Vin;
dist[12 * Np + n] = W2*Vin;
dist[13 * Np + n] = W2*Vin;
dist[16 * Np + n] = W2*Vin;
dist[17 * Np + n] = W2*Vin;
}
}
extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list,
double *dist,
double Vout,
int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
for (int idx = 0; idx < count; idx++) {
int n = list[idx];
dist[5 * Np + n] = W1*Vout;
dist[11 * Np + n] = W2*Vout;
dist[14 * Np + n] = W2*Vout;
dist[15 * Np + n] = W2*Vout;
dist[18 * Np + n] = W2*Vout;
}
}
extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList,
int *list,
double *dist,
double Vin, int count,
int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int nr5, nr11, nr14, nr15, nr18;
for (int idx = 0; idx < count; idx++) {
int n = list[idx];
// Unknown distributions
nr5 = d_neighborList[n + 4 * Np];
nr11 = d_neighborList[n + 10 * Np];
nr15 = d_neighborList[n + 14 * Np];
nr14 = d_neighborList[n + 13 * Np];
nr18 = d_neighborList[n + 17 * Np];
dist[nr5] = W1*Vin;
dist[nr11] = W2*Vin;
dist[nr15] = W2*Vin;
dist[nr14] = W2*Vin;
dist[nr18] = W2*Vin;
}
}
extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int nr6, nr12, nr13, nr16, nr17;
for (int idx = 0; idx < count; idx++) {
int n = list[idx];
// unknown distributions
nr6 = d_neighborList[n + 5 * Np];
nr12 = d_neighborList[n + 11 * Np];
nr16 = d_neighborList[n + 15 * Np];
nr17 = d_neighborList[n + 16 * Np];
nr13 = d_neighborList[n + 12 * Np];
dist[nr6] = W1*Vout;
dist[nr12] = W2*Vout;
dist[nr16] = W2*Vout;
dist[nr17] = W2*Vout;
dist[nr13] = W2*Vout;
}
}
extern "C" void ScaLBL_D3Q19_Poisson_Init(int *Map, double *dist, double *Psi,
int start, int finish, int Np) {
int n;

View File

@ -3,7 +3,7 @@
//#include <cuda_profiler_api.h>
#define NBLOCKS 1024
#define NTHREADS 256
#define NTHREADS 512
__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(int *neighborList,int *Map, double *dist, double *Psi, int start, int finish, int Np){
int n;
@ -328,7 +328,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map,
f16, f17, f18;
int nr1, nr2, nr3, nr4, nr5, nr6, nr7, nr8, nr9, nr10, nr11, nr12, nr13,
nr14, nr15, nr16, nr17, nr18;
double error,sum_q;
double sum_q;
double rlx = 1.0 / tau;
int idx;
@ -421,7 +421,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map,
f18 = dist[nr18];
sum_q = f1+f2+f3+f4+f5+f6+f7+f8+f9+f10+f11+f12+f13+f14+f15+f16+f17+f18;
error = 8.0*(sum_q - f0) + rho_e;
//error = 8.0*(sum_q - f0) + rho_e;
psi = 2.0*(f0*(1.0 - rlx) + rlx*(sum_q + 0.125*rho_e));
@ -545,6 +545,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist,
f17 = dist[18 * Np + n];
f18 = dist[17 * Np + n];
/* Ex = (f1 - f2) * rlx *
4.0; //NOTE the unit of electric field here is V/lu
Ey = (f3 - f4) * rlx *
4.0; //factor 4.0 is D3Q7 lattice squared speed of sound
Ez = (f5 - f6) * rlx * 4.0;
*/
Ex = (f1 - f2 + 0.5*(f7 - f8 + f9 - f10 + f11 - f12 + f13 - f14))*4.0; //NOTE the unit of electric field here is V/lu
Ey = (f3 - f4 + 0.5*(f7 - f8 - f9 + f10 + f15 - f16 + f17 - f18))*4.0;
Ez = (f5 - f6 + 0.5*(f11 - f12 - f13 + f14 + f15 - f16 - f17 + f18))*4.0;
@ -554,12 +560,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist,
sum_q = f1+f2+f3+f4+f5+f6+f7+f8+f9+f10+f11+f12+f13+f14+f15+f16+f17+f18;
error = 8.0*(sum_q - f0) + rho_e;
Error[n] = error;
psi = 2.0*(f0*(1.0 - rlx) + rlx*(sum_q + 0.125*rho_e));
idx = Map[n];
Psi[idx] = psi;
idx = Map[n];
Psi[idx] = psi;
// q = 0
dist[n] = W0*psi;//
@ -593,7 +601,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist,
dist[16 * Np + n] = W2*psi;//f16 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e;
dist[17 * Np + n] = W2*psi;//f17 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e;
dist[18 * Np + n] = W2*psi;//f18 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e;
//........................................................................
}
}
@ -635,6 +642,131 @@ __global__ void dvc_ScaLBL_D3Q19_Poisson_Init(int *Map, double *dist, double *P
}
}
__global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
int n = list[idx];
dist[6 * Np + n] = W1*Vin;
dist[12 * Np + n] = W2*Vin;
dist[13 * Np + n] = W2*Vin;
dist[16 * Np + n] = W2*Vin;
dist[17 * Np + n] = W2*Vin;
}
}
__global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
int n = list[idx];
dist[5 * Np + n] = W1*Vout;
dist[11 * Np + n] = W2*Vout;
dist[14 * Np + n] = W2*Vout;
dist[15 * Np + n] = W2*Vout;
dist[18 * Np + n] = W2*Vout;
}
}
__global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int nr5, nr11, nr14, nr15, nr18;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
int n = list[idx];
// Unknown distributions
nr5 = d_neighborList[n + 4 * Np];
nr11 = d_neighborList[n + 10 * Np];
nr15 = d_neighborList[n + 14 * Np];
nr14 = d_neighborList[n + 13 * Np];
nr18 = d_neighborList[n + 17 * Np];
dist[nr5] = W1*Vin;
dist[nr11] = W2*Vin;
dist[nr15] = W2*Vin;
dist[nr14] = W2*Vin;
dist[nr18] = W2*Vin;
}
}
__global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int nr6, nr12, nr13, nr16, nr17;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
int n = list[idx];
// unknown distributions
nr6 = d_neighborList[n + 5 * Np];
nr12 = d_neighborList[n + 11 * Np];
nr16 = d_neighborList[n + 15 * Np];
nr17 = d_neighborList[n + 16 * Np];
nr13 = d_neighborList[n + 12 * Np];
dist[nr6] = W1*Vout;
dist[nr12] = W2*Vout;
dist[nr16] = W2*Vout;
dist[nr17] = W2*Vout;
dist[nr13] = W2*Vout;
}
}
/* wrapper functions to launch kernels */
extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z<<<GRID,512>>>(list, dist, Vin, count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err));
}
}
//
extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np){
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z<<<GRID,512>>>(list, dist, Vout, count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count,int Np) {
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z<<<GRID,512>>>(d_neighborList, list, dist, Vin, count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err));
}
}
//
extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) {
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z<<<GRID,512>>>(d_neighborList, list, dist, Vout, count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map,
double *dist, double *Den_charge,
double *Psi, double *ElectricField,

View File

@ -0,0 +1,213 @@
********************************
Membrane Charging Dynamics
********************************
In this example, we consider membrane charging dynamics for a simple cell.
For the case considered in ``example/SingleCell`` an input membrane geometry is provided in the
file ``Bacterium.swc``, which specifies an oblong cell shape, relying on the ``.swc`` file format that
is commonly used to approximate neuron structures. The case considered is the four ion membrane transport
problem considered in Figure 4 from McClure & Li
The cell simulation is performed by the executable ``lbpm_nernst_planck_cell_simulator``, which is launched
in the same way as other parallel tools
.. code:: bash
mpirun -np 2 $LBPM_BIN/lbpm_nernst_planck_cell_simulator Bacterium.db
The input file ``Bacterium.db`` specifies the following
.. code:: c
MultiphysController {
timestepMax = 25000
num_iter_Ion_List = 4
analysis_interval = 100
tolerance = 1.0e-9
visualization_interval = 1000 // Frequency to write visualization data
}
Ions {
use_membrane = true
Restart = false
MembraneIonConcentrationList = 150.0e-3, 10.0e-3, 15.0e-3, 155.0e-3 //user-input unit: [mol/m^3]
temperature = 293.15 //unit [K]
number_ion_species = 4 //number of ions
tauList = 1.0, 1.0, 1.0, 1.0
IonDiffusivityList = 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec]
IonValenceList = 1, -1, 1, -1 //valence charge of ions; dimensionless; positive/negative integer
IonConcentrationList = 4.0e-3, 20.0e-3, 16.0e-3, 0.0e-3 //user-input unit: [mol/m^3]
BC_Solid = 0 //solid boundary condition; 0=non-flux BC; 1=surface ion concentration
FluidVelDummy = 0.0, 0.0, 0.0 // dummy fluid velocity for debugging
BC_InletList = 0, 0, 0, 0
BC_OutletList = 0, 0, 0, 0
}
Poisson {
lattice_scheme = "D3Q19"
epsilonR = 78.5 //fluid dielectric constant [dimensionless]
BC_Inlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
BC_Outlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
//--------------------------------------------------------------------------
//--------------------------------------------------------------------------
BC_Solid = 2 //solid boundary condition; 1=surface potential; 2=surface charge density
SolidLabels = 0 //solid labels for assigning solid boundary condition
SolidValues = 0 //if surface potential, unit=[V]; if surface charge density, unit=[C/m^2]
WriteLog = true //write convergence log for LB-Poisson solver
//------------------------------ advanced setting ------------------------------------
timestepMax = 4000 //max timestep for obtaining steady-state electrical potential
analysis_interval = 25 //timestep checking steady-state convergence
tolerance = 1.0e-10 //stopping criterion for steady-state solution
InitialValueLabels = 1, 2
InitialValues = 0.0, 0.0
}
Domain {
Filename = "Bacterium.swc"
nproc = 2, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 64, 64, 64 // Size of local domain (Nx,Ny,Nz)
N = 128, 64, 64 // size of the input image
voxel_length = 0.01 //resolution; user-input unit: [um]
BC = 0 // Boundary condition type
ReadType = "swc"
ReadValues = 0, 1, 2
WriteValues = 0, 1, 2
}
Analysis {
analysis_interval = 100
subphase_analysis_interval = 50 // Frequency to perform analysis
restart_interval = 5000 // Frequency to write restart data
restart_file = "Restart" // Filename to use for restart file (will append rank)
N_threads = 4 // Number of threads to use
load_balance = "independent" // Load balance method to use: "none", "default", "independent"
}
Visualization {
save_electric_potential = true
save_concentration = true
save_velocity = false
}
Membrane {
MembraneLabels = 2
VoltageThreshold = 0.0, 0.0, 0.0, 0.0
MassFractionIn = 1e-1, 1.0, 5e-3, 0.0
MassFractionOut = 1e-1, 1.0, 5e-3, 0.0
ThresholdMassFractionIn = 1e-1, 1.0, 5e-3, 0.0
ThresholdMassFractionOut = 1e-1, 1.0, 5e-3, 0.0
}
*******************
Example Output
*******************
Successful output looks like the following
.. code:: bash
********************************************************
Running LBPM Nernst-Planck Membrane solver
********************************************************
.... Read membrane permeability (MassFractionIn)
.... Read membrane permeability (MassFractionOut)
.... Read membrane permeability (ThresholdMassFractionIn)
.... Read membrane permeability (ThresholdMassFractionOut)
.... Read MembraneIonConcentrationList
voxel length = 0.010000 micron
voxel length = 0.010000 micron
Reading SWC file...
Number of lines in SWC file: 7
Number of lines extracted is: 7
shift swc data by 0.150000, 0.140000, 0.140000
Media porosity = 1.000000
LB Ion Solver: Initialized solid phase & converting to Signed Distance function
Domain set.
LB Ion Solver: Create ScaLBL_Communicator
LB Ion Solver: Set up memory efficient layout
LB Ion Solver: Allocating distributions
LB Ion Solver: Setting up device map and neighbor list
**** Creating membrane data structure ******
Number of active lattice sites (rank = 0): 262160
Membrane labels: 1
label=2, volume fraction = 0.133917
Creating membrane data structure...
Copy initial neighborlist...
Cut membrane links...
(cut 7105 links crossing membrane)
Construct membrane data structures...
Create device data structures...
Construct communication data structures...
Ion model setup complete
Analyze system with sub-domain size = 66 x 66 x 66
Set up analysis routines for 4 ions
LB Ion Solver: initializing D3Q7 distributions
...initializing based on membrane list
.... Set concentration(0): inside=0.15 [mol/m^3], outside=0.004 [mol/m^3]
.... Set concentration(1): inside=0.01 [mol/m^3], outside=0.02 [mol/m^3]
.... Set concentration(2): inside=0.015 [mol/m^3], outside=0.016 [mol/m^3]
.... Set concentration(3): inside=0.155 [mol/m^3], outside=0 [mol/m^3]
LB Ion Solver: initializing charge density
LB Ion Solver: solid boundary: non-flux boundary is assigned
LB Ion Solver: inlet boundary for Ion 1 is periodic
LB Ion Solver: outlet boundary for Ion 1 is periodic
LB Ion Solver: inlet boundary for Ion 2 is periodic
LB Ion Solver: outlet boundary for Ion 2 is periodic
LB Ion Solver: inlet boundary for Ion 3 is periodic
LB Ion Solver: outlet boundary for Ion 3 is periodic
LB Ion Solver: inlet boundary for Ion 4 is periodic
LB Ion Solver: outlet boundary for Ion 4 is periodic
*****************************************************
LB Ion Transport Solver:
Ion 1: LB relaxation tau = 1
Time conversion factor: 1.25e-08 [sec/lt]
Internal iteration: 2 [lt]
Ion 2: LB relaxation tau = 1
Time conversion factor: 1.25e-08 [sec/lt]
Internal iteration: 2 [lt]
Ion 3: LB relaxation tau = 1
Time conversion factor: 1.25e-08 [sec/lt]
Internal iteration: 2 [lt]
Ion 4: LB relaxation tau = 1
Time conversion factor: 1.25e-08 [sec/lt]
Internal iteration: 2 [lt]
*****************************************************
Ion model initialized
Main loop time_conv computed from ion 1: 2.5e-08[s/lt]
Main loop time_conv computed from ion 2: 2.5e-08[s/lt]
Main loop time_conv computed from ion 3: 2.5e-08[s/lt]
Main loop time_conv computed from ion 4: 2.5e-08[s/lt]
***********************************************************************************
LB-Poisson Solver: steady-state MaxTimeStep = 4000; steady-state tolerance = 1e-10
LB relaxation tau = 3.5
***********************************************************************************
LB-Poisson Solver: Use averaged MSE to check solution convergence.
LB-Poisson Solver: Use D3Q19 lattice structure.
voxel length = 0.010000 micron
voxel length = 0.010000 micron
Reading SWC file...
Number of lines in SWC file: 7
Number of lines extracted is: 7
shift swc data by 0.150000, 0.140000, 0.140000
Media porosity = 1.000000
LB-Poisson Solver: Initialized solid phase & converting to Signed Distance function
Domain set.
LB-Poisson Solver: Create ScaLBL_Communicator
LB-Poisson Solver: Set up memory efficient layout
LB-Poisson Solver: Allocating distributions
LB-Poisson Solver: Setting up device map and neighbor list
.... LB-Poisson Solver: check neighbor list
.... LB-Poisson Solver: copy neighbor list to GPU
Poisson solver created
LB-Poisson Solver: initializing D3Q19 distributions
LB-Poisson Solver: number of Poisson solid labels: 1
label=0, surface potential=0 [V], volume fraction=0
LB-Poisson Solver: number of Poisson initial-value labels: 2
label=1, initial potential=0 [V], volume fraction=0.96
label=2, initial potential=0 [V], volume fraction=0.13
POISSON MODEL: Reading restart file!
Poisson solver initialized
... getting Poisson solver error
-------------------------------------------------------------------
set coefficients
********************************************************
CPU time = 0.008526
Lattice update rate (per core)= 30.749833 MLUPS
Lattice update rate (total)= 61.499666 MLUPS
********************************************************

View File

@ -20,3 +20,5 @@ a basic introduction to working with LBPM.
color/*
analysis/*
membrane/*

View File

@ -1,4 +1,4 @@
\============
============
Publications
============
@ -13,3 +13,5 @@ Publications
* J.E. McClure, M. Berrill, W. Gray, C.T. Miller, C.T. "Tracking interface and common curve dynamics for two-fluid flow in porous media. Journal of Fluid Mechanics, 796, 211-232 (2016) https://doi.org/10.1017/jfm.2016.212
* J.E.McClure, D.Adalsteinsson, C.Pan, W.G.Gray, C.T.Miller "Approximation of interfacial properties in multiphase porous medium systems" Advances in Water Resources, 30 (3): 354-365 (2007) https://doi.org/10.1016/j.advwatres.2006.06.010
* J.E. McClure, Z. Li "Capturing membrane structure and function in lattice Boltzmann models" arXiv preprint arXiv:2208.14122 https://arxiv.org/pdf/2208.14122.pdf

View File

@ -1,6 +0,0 @@
=============================================
Poisson-Boltzmann model
=============================================
The LBPM Poisson-Boltzmann solver is designed to solve the Poisson-Boltzmann equation
to solve for the electric field in an ionic fluid.

View File

@ -0,0 +1,366 @@
=============================================
Cell model
=============================================
LBPM includes a whole-cell simulator based on a coupled solution of the Nernst-Planck equations with Gauss's law.
The resulting model is fully non-equilibrium, and can resolve the dynamics of how ions diffuse through the cellular
environment when subjected to complex membrane responses.
The lattice Boltzmann formulation is described below.
*********************
Nernst-Planck model
*********************
The Nernst-Planck model is designed to model ion transport based on the
Nernst-Planck equation.
.. math::
:nowrap:
$$
\frac{\partial C_k}{\partial t} + \nabla \cdot \mathbf{j}_k = 0
$$
where
.. math::
:nowrap:
$$
\mathbf{j}_k = C_k \mathbf{u} - D_k \Big( \nabla C_k + \frac{z_k C_k}{V_T} \nabla \psi\Big)
$$
A LBM solution is developed using a three-dimensional, seven velocity (D3Q7) lattice structure for each species. Each distribution is associated with a particular discrete velocity, such that the concentration is given by their sum,
.. math::
:nowrap:
$$
C_k = \sum_{q=0}^{6} f^k_q \;.
$$
Lattice Boltzmann equations (LBEs) are defined to determine the evolution of the distributions
.. math::
:nowrap:
$$
f^{k}_q (\mathbf{x}_n + \bm{\xi}_q \Delta t, t+ \Delta t)-
f^{k}_q (\mathbf{x}_n, t) = \frac{1}{\lambda_k}
\Big( f^{k}_q - f^{eq}_q \Big)\;,
$$
where the relaxation time :math:`\lambda_k` controls the bulk diffusion coefficient,
.. math::
:nowrap:
$$
D_k = c_s^2\Big( \lambda_k - \frac 12\Big)\;.
$$
The speed of sound for the D3Q7 lattice model is :math:`c_s^2 = \frac 14` and the weights are :math:`W_0 = 1/4` and :math:`W_1,\ldots, W_6 = 1/8`.
Equilibrium distributions are established from the fact that molecular velocity distribution follows a Gaussian distribution within the bulk fluids,
.. math::
:nowrap:
$$
f^{eq}_q = W_q C_k \Big[ 1 + \frac{\bm{\xi_q}\cdot \mathbf{u}^\prime}{c_s^2} \Big]\;.
$$
The velocity is given by
.. math::
:nowrap:
$$
\mathbf{u}^\prime = \mathbf{u} - \frac{z_k D_k}{V_T} \nabla \psi \;.
$$
Keys for the Nernst-Planck solver are provided in the ``Ion`` section of the input file database. Supported keys are
- ``use_membrane`` -- set up a membrane structure (defaults to ``true`` if not specified)
- ``Restart`` -- read concentrations from restart file (defaults to ``false`` if not specified)
- ``number_ion_species`` -- number of ions to use in the model
- ``temperature`` -- temperature to use for the thermal voltage (:math:`V_T=k_B T / e`, where the electron charge is :math:`e=1.6\times10^{-19}` Coulomb)
- ``FluidVelDummy`` -- vector providing a dummy fluid velocity field (for advection component)
- ``ElectricFieldDummy`` -- vectory providing a dummy electric field (for force component)
- ``tauList`` -- list of relaxation times to set the diffusion coefficient based on :math:`\lambda_k`.
- ``IonDiffusivityList`` -- list of physical ion diffusivities in units :math:`\mbox{m}^2/\mbox{second}`.
- ``IonValenceList`` -- list of ion valence charges for each ion in the model.
- ``IonConcentrationList`` -- list of concentrations to set for each ion.
- ``MembraneIonConcentrationList`` -- list of concentrations to set for each ion inside the membrane.
- ``BC_InletList`` -- boundary conditions for each ion at the z-inlet (``0`` for periodic, ``1`` to set concentration)
- ``BC_OutletList`` -- boundary conditions for each ion at the z-outlet
- ``InletValueList`` -- concentration value to set at the inlet (if not periodic)
- ``OutletValueList`` -- concentration value to set at the outlet (if not periodic)
*********************
Gauss's Law Model
*********************
The LBPM Gauss's law solver is designed to solve for the electric field in an ionic fluid.
.. math::
:nowrap:
$$
\nabla^2_{fe} \psi (\mathbf{x}_i) = \frac{1}{6 \Delta x^2}
\Bigg( 2 \sum_{q=1}^{6} \psi(\mathbf{x}_i + \bm{\xi}_q \Delta t)
+ \sum_{q=7}^{18} \psi(\mathbf{x}_i + \bm{\xi}_q \Delta t)
- 24 \psi (\mathbf{x}_i) \Bigg) \;,
$$
The equilibrium functions are defined as
.. math::
:nowrap:
$$
g_q^{eq} = w_q \psi\;,
$$
where :math:`w_0=1/2`, :math:`w_q=1/24` for :math:`q=1,\ldots,6` and :math:`w_q=1/48` for :math:`q=7,\ldots,18`
which implies that
.. math::
:nowrap:
$$
\psi = \sum_{q=0}^{Q} g_q^{eq}\;.
$$
Given a particular initial condition for :math:`\psi`, let us consider application of the standard D3Q19 streaming step based on the equilibrium distributions
.. math::
:nowrap:
$$
g_q^\prime(\mathbf{x}, t) = g_q^{eq}(\mathbf{x}-\bm{\xi}_q\Delta t, t+ \Delta t)\;.
$$
Relative to the solution of Gauss's law, the error is given by
.. math::
:nowrap:
$$
\varepsilon_{\psi} =
8 \Big[ -g_0 + \sum_{q=1}^Q g_q^\prime(\mathbf{x}, t) \Big]
+ \frac{\rho_e}{\epsilon_r \epsilon_0} \;.
$$
Using the fact that :math:`f_0 = W_0 \psi`, we can compute the value
:math:`\psi^\prime` that would kill the error. We set :math:`\varepsilon_{\psi}=0`
and rearrange terms to obtain
.. math::
:nowrap:
$$
\psi^\prime (\mathbf{x},t) = \frac{1}{W_0}\Big[ \sum_{q=1}^Q g_q^\prime(\mathbf{x}, t)
+ \frac{1}{8}\frac{\rho_e}{\epsilon_r \epsilon_0}\Big] \;.
$$
The local value of the potential is then updated based on a relaxation scheme, which is controlled by the relaxation time :math:`\tau_\psi`
.. math::
:nowrap:
$$
\psi(\mathbf{x},t+\Delta t) \leftarrow \Big(1 - \frac{1}{\tau_\psi} \Big )\psi (\mathbf{x},t)
+ \frac{1}{\tau_\psi} \psi^\prime (\mathbf{x},t)\;.
$$
The algorithm can then proceed to the next timestep.
Keys to control the Gauss's law solver are specified in the ``Poisson`` section of the input database.
Supported keys are:
- ``Restart`` -- read electric potential from a restart file (default ``false``)
- ``timestepMax`` -- maximum number of timesteps to run before exiting
- ``tau`` -- relaxation time
- ``analysis_interval`` -- how often to check solution for steady state
- ``tolerance`` -- controls the required accuracy
- ``epsilonR`` -- controls the electric permittivity
- ``WriteLog`` -- write a convergence log
***************************
Membrane Model
***************************
The LBPM membrane model provides the basis to model cellular dynamics.
There are currently two supported ways to specify the membrane location:
1. provide a segemented image that is labeled to differentiate the cell
interior and exterior. See the script ``NaCl-cell.py`` and input file ``NaCl.db`` as a reference for how to use labeled images.
- ``IonConcentrationFile`` -- list of files that specify the initial concentration for each ion
- ``Filename`` -- 8-bit binary file provided in the ``Domain`` section of the input database
- ``ReadType`` -- this should be ``"8bit"`` (this is the default)
2. provide a ``.swc`` file that specifies the geometry (see example input file below).
- ``Filename`` -- swc file name should be provided in the ``Domain`` section of the input database
- ``ReadType`` -- this should be ``"swc"`` (required since ``"8bit"`` is the internal default)
Example input files for both cases are stored within the LBPM repository, located at ``example/SingleCell/``
The membrane simply prevents the diffusion of ions. All lattice links crossing the membrane are stored in a dedicated data structure so that transport is decoupled from the bulk regions. Suppose that site :math:`\mathbf{x}_{q\ell}` is inside the membrane and :math:`\mathbf{x}_{p\ell}` is outside the membrane, with :math:`\mathbf{x}_{p \ell } = \mathbf{x}_{q\ell} + \bm{\xi}_q \Delta t`. For each species :math:`k`, transport across each link :math:`\ell` is controlled by a pair of coefficients, :math:`\alpha^k_{\ell p}` and :math:`\alpha^k_{\ell q}`. Ions transported from the outside to the inside are transported by the particular distribution that is associated with the direction :math:`\xi_q`
.. math::
:nowrap:
$$
{ f_{q}^{k \prime} (\mathbf{x}_{q \ell}) \gets (1-\alpha^k_{\ell q}) f_{q}^{k} (\mathbf{x}_{q\ell}) + \alpha^k_{\ell p } f_{ p}^{k} (\mathbf{x}_{p\ell})}
$$
Similarly, for ions transported from the inside to the outside
.. math::
:nowrap:
$$
{f_{p}^{k \prime} (\mathbf{x}_{p\ell}) \gets (1-\alpha^k_{\ell p}) f_{p}^{k} (\mathbf{x}_{p\ell}) + \alpha^k_{\ell q } f_{q}^{k} (\mathbf{x}_{q\ell})}
$$
The basic closure relationship that is implemented is for voltage-gated ion channels.
Let :math:`\Delta \psi_\ell = \psi(\mathbf{x}_{p\ell} ,t) - \psi(\mathbf{x}_{q\ell},t)` be the membrane potential across link :math:`\ell`. Since :math:`\psi` is determined based on the charge density, :math:`\Delta \psi_\ell` can vary with both space and time. The behavior of the gate is implmented as follows,
.. math::
:nowrap:
$$
\Delta \psi_\ell > \tilde{V}_m\; \Rightarrow \; \mbox{gate is open} \; \Rightarrow \; \alpha^{k}_{q \ell} = \alpha_{1} + \alpha_2\;,
$$
and
.. math::
:nowrap:
$$
\Delta \psi_\ell \le \tilde{V}_m\; \Rightarrow \; \mbox{gate is closed}\; \Rightarrow \; \alpha^{{k}}_{q \ell} = \alpha_1\;
$$
where :math:`\tilde{V}_m` is the membrane voltage threshold that controls gate. Mass conservation dictates that
.. math::
:nowrap:
$$
\alpha_1 \ge 0\;, \quad \alpha_2 \ge 0\;, \quad \alpha_1 + \alpha_2 \le 1\;.
$$
The rule is enforced based on the Heaviside function, as follows
.. math::
:nowrap:
$$
\alpha_{\ell q}^{k} (\Delta \psi_\ell) = \alpha_1 + \alpha_2 H\big(\Delta \psi_\ell - \tilde{V}_m \big)\;.
$$
Note that different coefficients are specified for each ion in the model.
Keys for the membrane model are set in the ``Membrane`` section of the input file database. Supported keys are
- ``VoltageThreshold`` -- voltage threshold (may be different for each ion)
- ``MassFractionIn`` -- value of :math:`\alpha^k_{\ell p}` when the voltage threshold is not met
- ``MassFractionOut`` -- value of :math:`\alpha^k_{\ell q}` when the voltage threshold is not met
- ``ThresholdMassFractionIn`` -- value of :math:`\alpha^k_{\ell p}` when the voltage threshold is met
- ``ThresholdMassFractionOut`` -- value of :math:`\alpha^k_{\ell q}` when the voltage threshold is met
****************************
Example Input File
****************************
.. code-block:: c
MultiphysController {
timestepMax = 25000
num_iter_Ion_List = 4
analysis_interval = 100
tolerance = 1.0e-9
visualization_interval = 1000 // Frequency to write visualization data
}
Ions {
use_membrane = true
Restart = false
MembraneIonConcentrationList = 150.0e-3, 10.0e-3, 15.0e-3, 155.0e-3 //user-input unit: [mol/m^3]
temperature = 293.15 //unit [K]
number_ion_species = 4 //number of ions
tauList = 1.0, 1.0, 1.0, 1.0
IonDiffusivityList = 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec]
IonValenceList = 1, -1, 1, -1 //valence charge of ions; dimensionless; positive/negative integer
IonConcentrationList = 4.0e-3, 20.0e-3, 16.0e-3, 0.0e-3 //user-input unit: [mol/m^3]
BC_Solid = 0 //solid boundary condition; 0=non-flux BC; 1=surface ion concentration
//SolidLabels = 0 //solid labels for assigning solid boundary condition; ONLY for BC_Solid=1
//SolidValues = 1.0e-5 // user-input surface ion concentration unit: [mol/m^2]; ONLY for BC_Solid=1
FluidVelDummy = 0.0, 0.0, 0.0 // dummy fluid velocity for debugging
BC_InletList = 0, 0, 0, 0
BC_OutletList = 0, 0, 0, 0
}
Poisson {
lattice_scheme = "D3Q19"
epsilonR = 78.5 //fluid dielectric constant [dimensionless]
BC_Inlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
BC_Outlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
//--------------------------------------------------------------------------
//--------------------------------------------------------------------------
BC_Solid = 2 //solid boundary condition; 1=surface potential; 2=surface charge density
SolidLabels = 0 //solid labels for assigning solid boundary condition
SolidValues = 0 //if surface potential, unit=[V]; if surface charge density, unit=[C/m^2]
WriteLog = true //write convergence log for LB-Poisson solver
// ------------------------------- Testing Utilities ----------------------------------------
// ONLY for code debugging; the followings test sine/cosine voltage BCs; disabled by default
TestPeriodic = false
TestPeriodicTime = 1.0 //unit:[sec]
TestPeriodicTimeConv = 0.01 //unit:[sec]
TestPeriodicSaveInterval = 0.2 //unit:[sec]
//------------------------------ advanced setting ------------------------------------
timestepMax = 4000 //max timestep for obtaining steady-state electrical potential
analysis_interval = 25 //timestep checking steady-state convergence
tolerance = 1.0e-10 //stopping criterion for steady-state solution
InitialValueLabels = 1, 2
InitialValues = 0.0, 0.0
}
Domain {
Filename = "Bacterium.swc"
nproc = 2, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 64, 64, 64 // Size of local domain (Nx,Ny,Nz)
N = 128, 64, 64 // size of the input image
voxel_length = 0.01 //resolution; user-input unit: [um]
BC = 0 // Boundary condition type
ReadType = "swc"
ReadValues = 0, 1, 2
WriteValues = 0, 1, 2
}
Analysis {
analysis_interval = 100
subphase_analysis_interval = 50 // Frequency to perform analysis
restart_interval = 5000 // Frequency to write restart data
restart_file = "Restart" // Filename to use for restart file (will append rank)
N_threads = 4 // Number of threads to use
load_balance = "independent" // Load balance method to use: "none", "default", "independent"
}
Visualization {
save_electric_potential = true
save_concentration = true
save_velocity = false
}
Membrane {
MembraneLabels = 2
VoltageThreshold = 0.0, 0.0, 0.0, 0.0
MassFractionIn = 1e-1, 1.0, 5e-3, 0.0
MassFractionOut = 1e-1, 1.0, 5e-3, 0.0
ThresholdMassFractionIn = 1e-1, 1.0, 5e-3, 0.0
ThresholdMassFractionOut = 1e-1, 1.0, 5e-3, 0.0
}

View File

@ -176,42 +176,70 @@ The non-zero equilibrium moments are defined as
:nowrap:
$$
m_1^{eq} = (j_x^2+j_y^2+j_z^2) - \alpha |\textbf{C}|, \\
m_1^{eq} = 19\frac{ j_x^2+j_y^2+j_z^2}{\rho_0} - 11\rho - 19 \alpha |\textbf{C}|, \\
$$
.. math::
:nowrap:
$$
m_2^{eq} = 3\rho - \frac{11( j_x^2+j_y^2+j_z^2)}{2\rho_0}, \\
$$
.. math::
:nowrap:
$$
m_4^{eq} = -\frac{2 j_x}{3}, \\
$$
.. math::
:nowrap:
$$
m_6^{eq} = -\frac{2 j_y}{3}, \\
$$
.. math::
:nowrap:
$$
m_8^{eq} = -\frac{2 j_z}{3}, \\
$$
.. math::
:nowrap:
$$
m_9^{eq} = (2j_x^2-j_y^2-j_z^2)+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\
m_9^{eq} = \frac{2j_x^2-j_y^2-j_z^2}{\rho_0}+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{11}^{eq} = (j_y^2-j_z^2) + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\
m_{11}^{eq} = \frac{j_y^2-j_z^2}{\rho_0} + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{13}^{eq} = j_x j_y + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\
m_{13}^{eq} = \frac{j_x j_y}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\
$$
.. math::
:nowrap:
$$
m_{14}^{eq} = j_y j_z + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\
m_{14}^{eq} = \frac{j_y j_z}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\
$$
.. math::
:nowrap:
$$
m_{15}^{eq} = j_x j_z + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;,
m_{15}^{eq} = \frac{j_x j_z}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;.
$$
where the color gradient is determined from the phase indicator field

View File

@ -241,46 +241,75 @@ The relaxation parameters are determined from the relaxation time:
The non-zero equilibrium moments are defined as
.. math::
:nowrap:
$$
m_1^{eq} = (j_x^2+j_y^2+j_z^2) - \alpha |\textbf{C}|, \\
m_1^{eq} = 19\frac{ j_x^2+j_y^2+j_z^2}{\rho_0} - 11\rho - 19 \alpha |\textbf{C}|, \\
$$
.. math::
:nowrap:
$$
m_2^{eq} = 3\rho - \frac{11( j_x^2+j_y^2+j_z^2)}{2\rho_0}, \\
$$
.. math::
:nowrap:
$$
m_4^{eq} = -\frac{2 j_x}{3}, \\
$$
.. math::
:nowrap:
$$
m_6^{eq} = -\frac{2 j_y}{3}, \\
$$
.. math::
:nowrap:
$$
m_8^{eq} = -\frac{2 j_z}{3}, \\
$$
.. math::
:nowrap:
$$
m_9^{eq} = (2j_x^2-j_y^2-j_z^2)+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\
m_9^{eq} = \frac{2j_x^2-j_y^2-j_z^2}{\rho_0}+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{11}^{eq} = (j_y^2-j_z^2) + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\
m_{11}^{eq} = \frac{j_y^2-j_z^2}{\rho_0} + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\
$$
.. math::
:nowrap:
$$
m_{13}^{eq} = j_x j_y + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\
m_{13}^{eq} = \frac{j_x j_y}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\
$$
.. math::
:nowrap:
$$
m_{14}^{eq} = j_y j_z + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\
m_{14}^{eq} = \frac{j_y j_z}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\
$$
.. math::
:nowrap:
$$
m_{15}^{eq} = j_x j_z + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;,
m_{15}^{eq} = \frac{j_x j_z}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;.
$$
where the color gradient is determined from the phase indicator field

View File

@ -12,9 +12,7 @@ Currently supported lattice Boltzmann models
mrt/*
nernstPlanck/*
PoissonBoltzmann/*
cell/*
greyscale/*

View File

@ -1,6 +0,0 @@
=============================================
Nernst-Planck model
=============================================
The Nernst-Planck model is designed to model ion transport based on the
Nernst-Planck equation.

View File

@ -1,41 +1,38 @@
import numpy as np
import matplotlib.pylab as plt
Nx = 64
Ny = 64
Nz = 64
cx = Nx/2
cy = Ny/2
cz = Nz/2
radius = 12
D=np.ones((40,40,40),dtype="uint8")
D=np.ones((Nx,Ny,Nz),dtype="uint8")
cx = 20
cy = 20
cz = 20
radius = 8
for i in range(0, Nx):
for j in range (0, Ny):
for k in range (0,Nz):
for i in range(0, 40):
for j in range (0, 40):
for k in range (0,40):
dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
if (dist < radius ) :
D[i,j,k] = 2
D.tofile("cell_64x64x64.raw")
D.tofile("cell_40x40x40.raw")
C1=np.zeros((Nx,Ny,Nz),dtype="double")
C2=np.zeros((Nx,Ny,Nz),dtype="double")
C1=np.zeros((40,40,40),dtype="double")
C2=np.zeros((40,40,40),dtype="double")
for i in range(0, Nx):
for j in range (0, Ny):
for k in range (0,Nz):
for i in range(0, 40):
for j in range (0, 40):
for k in range (0,40):
#outside the cell
C1[i,j,k] = 125.0e-6 # Na
C2[i,j,k] = 125.0e-6 # Cl
dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
# inside the cell
if (dist < radius ) :
C1[i,j,k] = 110.0e-6
C2[i,j,k] = 110.0e-6
C1[i,j,k] = 5.0e-6
C2[i,j,k] = 5.0e-6
C1.tofile("cell_concentration_Na_64x64x64.raw")
C2.tofile("cell_concentration_Cl_64x64x64.raw")
C1.tofile("cell_concentration_Na_40x40x40.raw")
C2.tofile("cell_concentration_Cl_40x40x40.raw")

View File

@ -1,9 +1,10 @@
MultiphysController {
timestepMax = 2000
timestepMax = 20
num_iter_Ion_List = 2
analysis_interval = 40
analysis_interval = 50
tolerance = 1.0e-9
visualization_interval = 40 // Frequency to write visualization data
visualization_interval = 100 // Frequency to write visualization data
analysis_interval = 50 // Frequency to perform analysis
}
Stokes {
tau = 1.0
@ -12,7 +13,9 @@ Stokes {
nu_phys = 0.889e-6 //fluid kinematic viscosity; user-input unit: [m^2/sec]
}
Ions {
IonConcentrationFile = "cell_concentration_Na_64x64x64.raw", "double", "cell_concentration_Cl_64x64x64.raw", "double"
use_membrane = true
Restart = false
IonConcentrationFile = "cell_concentration_Na_40x40x40.raw", "double", "cell_concentration_Cl_40x40x40.raw", "double"
temperature = 293.15 //unit [K]
number_ion_species = 2 //number of ions
tauList = 1.0, 1.0
@ -25,6 +28,8 @@ Ions {
FluidVelDummy = 0.0, 0.0, 0.0 // dummy fluid velocity for debugging
}
Poisson {
lattice_scheme = "D3Q19"
Restart = false
epsilonR = 78.5 //fluid dielectric constant [dimensionless]
BC_Inlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
BC_Outlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
@ -41,16 +46,19 @@ Poisson {
TestPeriodicTimeConv = 0.01 //unit:[sec]
TestPeriodicSaveInterval = 0.2 //unit:[sec]
//------------------------------ advanced setting ------------------------------------
timestepMax = 4000 //max timestep for obtaining steady-state electrical potential
analysis_interval = 25 //timestep checking steady-state convergence
tolerance = 1.0e-10 //stopping criterion for steady-state solution
timestepMax = 100000 //max timestep for obtaining steady-state electrical potential
analysis_interval = 200 //timestep checking steady-state convergence
tolerance = 1.0e-6 //stopping criterion for steady-state solution
InitialValueLabels = 1, 2
InitialValues = 0.0, 0.0
}
Domain {
Filename = "cell_64x64x64.raw"
Filename = "cell_40x40x40.raw"
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 64, 64, 64 // Size of local domain (Nx,Ny,Nz)
N = 64, 64, 64 // size of the input image
voxel_length = 0.1 //resolution; user-input unit: [um]
n = 40, 40, 40 // Size of local domain (Nx,Ny,Nz)
N = 40, 40, 40 // size of the input image
voxel_length = 1.0 //resolution; user-input unit: [um]
BC = 0 // Boundary condition type
ReadType = "8bit"
ReadValues = 0, 1, 2
@ -72,9 +80,9 @@ Visualization {
Membrane {
MembraneLabels = 2
VoltageThreshold = 0.0, 0.0
MassFractionIn = 1e-2, 1e-8
MassFractionOut = 1e-2, 1e-8
ThresholdMassFractionIn = 1e-2, 1e-8
ThresholdMassFractionOut = 1e-2, 1e-8
MassFractionIn = 0e-2, 5e-2
MassFractionOut = 0e-2, 5e-2
ThresholdMassFractionIn = 0e-2, 5e-2
ThresholdMassFractionOut = 0e-2, 5e-2
}

View File

@ -301,7 +301,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map,
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if (n<finish) {
//Load data
//Load data
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
@ -451,7 +451,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map,
dist[nr18] = W2*psi; //f17 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e;
// q = 18
dist[nr17] = W2*psi;
dist[nr17] = W2*psi; //f18 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e;
}
}
}
@ -479,7 +480,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist,
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if (n<finish) {
//Load data
//Load data
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
//and thus the net space charge density is zero.
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
@ -505,6 +506,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist,
f17 = dist[18 * Np + n];
f18 = dist[17 * Np + n];
/* Ex = (f1 - f2) * rlx *
4.0; //NOTE the unit of electric field here is V/lu
Ey = (f3 - f4) * rlx *
4.0; //factor 4.0 is D3Q7 lattice squared speed of sound
Ez = (f5 - f6) * rlx * 4.0;
*/
Ex = (f1 - f2 + 0.5*(f7 - f8 + f9 - f10 + f11 - f12 + f13 - f14))*4.0; //NOTE the unit of electric field here is V/lu
Ey = (f3 - f4 + 0.5*(f7 - f8 - f9 + f10 + f15 - f16 + f17 - f18))*4.0;
Ez = (f5 - f6 + 0.5*(f11 - f12 - f13 + f14 + f15 - f16 - f17 + f18))*4.0;
@ -514,12 +521,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist,
sum_q = f1+f2+f3+f4+f5+f6+f7+f8+f9+f10+f11+f12+f13+f14+f15+f16+f17+f18;
error = 8.0*(sum_q - f0) + rho_e;
Error[n] = error;
psi = 2.0*(f0*(1.0 - rlx) + rlx*(sum_q + 0.125*rho_e));
idx = Map[n];
Psi[idx] = psi;
idx = Map[n];
Psi[idx] = psi;
// q = 0
dist[n] = W0*psi;//
@ -553,7 +562,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist,
dist[16 * Np + n] = W2*psi;//f16 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e;
dist[17 * Np + n] = W2*psi;//f17 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e;
dist[18 * Np + n] = W2*psi;//f18 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e;
//........................................................................
}
}
@ -594,6 +602,129 @@ __global__ void dvc_ScaLBL_D3Q19_Poisson_Init(int *Map, double *dist, double *P
}
}
}
__global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
int n = list[idx];
dist[6 * Np + n] = W1*Vin;
dist[12 * Np + n] = W2*Vin;
dist[13 * Np + n] = W2*Vin;
dist[16 * Np + n] = W2*Vin;
dist[17 * Np + n] = W2*Vin;
}
}
__global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
int n = list[idx];
dist[5 * Np + n] = W1*Vout;
dist[11 * Np + n] = W2*Vout;
dist[14 * Np + n] = W2*Vout;
dist[15 * Np + n] = W2*Vout;
dist[18 * Np + n] = W2*Vout;
}
}
__global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int nr5, nr11, nr14, nr15, nr18;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
int n = list[idx];
// Unknown distributions
nr5 = d_neighborList[n + 4 * Np];
nr11 = d_neighborList[n + 10 * Np];
nr15 = d_neighborList[n + 14 * Np];
nr14 = d_neighborList[n + 13 * Np];
nr18 = d_neighborList[n + 17 * Np];
dist[nr5] = W1*Vin;
dist[nr11] = W2*Vin;
dist[nr15] = W2*Vin;
dist[nr14] = W2*Vin;
dist[nr18] = W2*Vin;
}
}
__global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) {
double W1 = 1.0/24.0;
double W2 = 1.0/48.0;
int nr6, nr12, nr13, nr16, nr17;
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx < count){
int n = list[idx];
// unknown distributions
nr6 = d_neighborList[n + 5 * Np];
nr12 = d_neighborList[n + 11 * Np];
nr16 = d_neighborList[n + 15 * Np];
nr17 = d_neighborList[n + 16 * Np];
nr13 = d_neighborList[n + 12 * Np];
dist[nr6] = W1*Vout;
dist[nr12] = W2*Vout;
dist[nr16] = W2*Vout;
dist[nr17] = W2*Vout;
dist[nr13] = W2*Vout;
}
}
/* wrapper functions to launch kernels */
extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z<<<GRID,512>>>(list, dist, Vin, count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("Hip error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err));
}
}
//
extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np){
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z<<<GRID,512>>>(list, dist, Vout, count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("Hip error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count,int Np) {
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z<<<GRID,512>>>(d_neighborList, list, dist, Vin, count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("Hip error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err));
}
}
//
extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) {
int GRID = count / 512 + 1;
dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z<<<GRID,512>>>(d_neighborList, list, dist, Vout, count, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("Hip error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map,
double *dist, double *Den_charge,

View File

@ -301,6 +301,7 @@ void ScaLBL_IonModel::ReadParams(string filename) {
//NOTE: the maximum iteration timesteps for ions are left unspecified
// it relies on the multiphys controller to compute the max timestep
USE_MEMBRANE = true;
Restart = false;
// read the input database
db = std::make_shared<Database>(filename);
domain_db = db->getDatabase("Domain");
@ -1080,20 +1081,42 @@ void ScaLBL_IonModel::Initialize() {
printf("LB Ion Solver: initializing D3Q7 distributions\n");
//USE_MEMBRANE = true;
if (USE_MEMBRANE){
double *Ci_host;
if (rank == 0)
printf(" ...initializing based on membrane list \n");
Ci_host = new double[number_ion_species * Np];
for (size_t ic = 0; ic < number_ion_species; ic++) {
AssignIonConcentrationMembrane( &Ci_host[ic * Np], ic);
}
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species * sizeof(double) * Np);
comm.barrier();
for (size_t ic = 0; ic < number_ion_species; ic++) {
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic * Np * 7], &Ci[ic * Np], Np);
}
delete[] Ci_host;
double *Ci_host;
Ci_host = new double[number_ion_species * Np];
if (ion_db->keyExists("IonConcentrationFile")) {
//NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
auto File_ion = ion_db->getVector<std::string>("IonConcentrationFile");
if (File_ion.size() == 2*number_ion_species) {
for (size_t ic = 0; ic < number_ion_species; ic++) {
AssignIonConcentration_FromFile(&Ci_host[ic * Np], File_ion,
ic);
}
ScaLBL_CopyToDevice(Ci, Ci_host,
number_ion_species * sizeof(double) * Np);
comm.barrier();
for (size_t ic = 0; ic < number_ion_species; ic++) {
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic * Np * 7], &Ci[ic * Np],
Np);
}
}
}
else{
if (rank == 0)
printf(" ...initializing based on membrane list \n");
for (size_t ic = 0; ic < number_ion_species; ic++) {
AssignIonConcentrationMembrane( &Ci_host[ic * Np], ic);
}
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species * sizeof(double) * Np);
comm.barrier();
for (size_t ic = 0; ic < number_ion_species; ic++) {
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic * Np * 7], &Ci[ic * Np], Np);
}
}
delete[] Ci_host;
}
else if (ion_db->keyExists("IonConcentrationFile")) {
//NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
auto File_ion = ion_db->getVector<std::string>("IonConcentrationFile");
@ -1549,7 +1572,6 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
IonMembrane->IonTransport(&fq[ic * Np * 7],&Ci[ic * Np]);
/* if (BoundaryConditionSolid == 1) {
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic * Np * 7], IonSolid);

View File

@ -392,19 +392,17 @@ void ScaLBL_MRTModel::Run() {
double h = Dm->voxel_length;
double absperm =
h * h * mu * Mask->Porosity() * flow_rate / force_mag;
double absperm_adjusted =
h * h * mu * Mask->Porosity() * Mask->Porosity() * flow_rate / force_mag;
absperm_adjusted *= 1013.0; // Convert to mDarcy
absperm *= 1013.0; // Convert to mDarcy
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 %.8g\n",
"%.8g %.8g\n",
timestep, Fx, Fy, Fz, mu, h * h * h * Vs, h * h * As,
h * Hs, Xs, vax, vay, vaz, absperm, absperm_adjusted);
h * Hs, Xs, vax, vay, vaz, absperm);
fclose(log_file);
}
}

View File

@ -86,7 +86,7 @@ void ScaLBL_Poisson::ReadParams(string filename){
}
//'tolerance_method' can be {"MSE","MSE_max"}
tolerance_method = electric_db->getWithDefault<std::string>( "tolerance_method", "MSE" );
lattice_scheme = electric_db->getWithDefault<std::string>( "lattice_scheme", "D3Q7" );
lattice_scheme = electric_db->getWithDefault<std::string>( "lattice_scheme", "D3Q19" );
if (electric_db->keyExists( "epsilonR" )){
epsilonR = electric_db->getScalar<double>( "epsilonR" );
}
@ -726,13 +726,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times
// *************ODD TIMESTEP*************//
timestep++;
SolveElectricPotentialAAodd(timestep_from_Study);//update electric potential
SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision
SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC,timestep);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier();
// *************EVEN TIMESTEP*************//
timestep++;
SolveElectricPotentialAAeven(timestep_from_Study);//update electric potential
SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision
SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC,timestep);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
@ -797,18 +797,17 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times
//************************************************************************/
// *************ODD TIMESTEP*************//
timestep++;
//SolveElectricPotentialAAodd(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential
SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision
//SolveElectricPotentialAAodd(timestep_from_Study);//,ChargeDensity, UseSlippingVelBC);//update electric potential
SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC,timestep);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier();
// *************EVEN TIMESTEP*************//
timestep++;
//SolveElectricPotentialAAeven(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential
SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision
//SolveElectricPotentialAAeven(timestep_from_Study);//,ChargeDensity, UseSlippingVelBC);//update electric potential
SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC,timestep);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
// Check convergence of steady-state solution
if (timestep==2){
//save electric potential for convergence check
@ -863,6 +862,144 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times
}
}
void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bool UseSlippingVelBC, int timestep_from_Study){
double error = 1.0;
double threshold = 10000000.0;
bool SET_THRESHOLD = false;
if (electric_db->keyExists( "rescale_at_distance" )){
SET_THRESHOLD = true;
threshold = electric_db->getScalar<double>( "rescale_at_distance" );
}
if (BoundaryConditionInlet > 0) SET_THRESHOLD = false;
if (BoundaryConditionOutlet > 0) SET_THRESHOLD = false;
double *host_Error;
host_Error = new double [Np];
timestep=0;
auto t1 = std::chrono::system_clock::now();
while (timestep < timestepMax && error > tolerance) {
//************************************************************************/
// *************ODD TIMESTEP*************//
timestep++;
//SolveElectricPotentialAAodd(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential
SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC,timestep);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier();
// *************EVEN TIMESTEP*************//
timestep++;
//SolveElectricPotentialAAeven(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential
SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC,timestep);//perform collision
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
// Check convergence of steady-state solution
if (timestep==2){
//save electric potential for convergence check
}
if (timestep%analysis_interval==0){
/* get the elecric potential */
ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
if (rank==0) printf(" ... getting Poisson solver error \n");
double err = 0.0;
double max_error = 0.0;
ScaLBL_CopyToHost(host_Error,ResidualError,sizeof(double)*Np);
for (int idx=0; idx<Np; idx++){
err = host_Error[idx]*host_Error[idx];
if (err > max_error ){
max_error = err;
}
}
error=Dm->Comm.maxReduce(max_error);
if (error > tolerance && SET_THRESHOLD){
/* don't use this with an external BC */
// cpompute the far-field electric potential
double inside_local = 0.0;
double outside_local = 0.0;
double inside_count_local = 0.0;
double outside_count_local = 0.0;
/* global values */
double inside_global = 0.0;
double outside_global = 0.0;
double inside_count_global = 0.0;
double outside_count_global = 0.0;
for (int k=1; k<Nz; k++){
for (int j=1; j<Ny; j++){
for (int i=1; i<Nx; i++){
int n = k*Nx*Ny + j*Nx + i;
double distance = MembraneDistance(i,j,k);
if (distance > threshold && distance < (threshold + 1.0)){
outside_count_local += 1.0;
outside_local += Psi_host(n);
}
else if (distance < (-1.0)*threshold && distance > (-1.0)*(threshold + 1.0)){
inside_count_local += 1.0;
inside_local += Psi_host(n);
}
}
}
}
inside_count_global = Dm->Comm.sumReduce(inside_count_local);
outside_count_global = Dm->Comm.sumReduce(outside_count_local);
outside_global = Dm->Comm.sumReduce(outside_local);
inside_global = Dm->Comm.sumReduce(inside_local);
outside_global /= outside_count_global;
inside_global /= inside_count_global;
if (rank==0) printf(" Rescaling far-field electric potential to value (outside): %f \n",outside_global);
if (rank==0) printf(" Rescaling far-field electric potential to value (inside): %f \n",inside_global);
// rescale the far-field electric potential
for (int k=1; k<Nz; k++){
for (int j=1; j<Ny; j++){
for (int i=1; i<Nx; i++){
int n = k*Nx*Ny + j*Nx + i;
double distance = MembraneDistance(i,j,k);
if ( distance > (threshold + 1.0)){
Psi_host(n) = outside_global;
}
else if ( distance < (-1.0)*(threshold + 1.0)){
Psi_host(n) = inside_global;
}
}
}
}
ScaLBL_CopyToDevice(Psi,Psi_host.data(),sizeof(double)*Nx*Ny*Nz);
}
/* compute the eletric field */
//ScaLBL_D3Q19_Poisson_getElectricField(fq, ElectricField, tau, Np);
}
}
if (rank == 0)
printf("---------------------------------------------------------------"
"----\n");
// Compute the walltime per timestep
auto t2 = std::chrono::system_clock::now();
double cputime = std::chrono::duration<double>(t2 - t1).count() / 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");
delete [] host_Error;
//************************************************************************/
if(WriteLog==true){
getConvergenceLog(timestep,error);
}
}
void ScaLBL_Poisson::getConvergenceLog(int timestep,double error){
if ( rank == 0 ) {
fprintf(TIMELOG,"%i %.5g\n",timestep,error);
@ -913,22 +1050,22 @@ void ScaLBL_Poisson::SolveElectricPotentialAAodd(int timestep_from_Study){
if (BoundaryConditionInlet > 0){
switch (BoundaryConditionInlet){
case 1:
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study);
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}
if (BoundaryConditionOutlet > 0){
switch (BoundaryConditionOutlet){
case 1:
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}
@ -983,22 +1120,22 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){
if (BoundaryConditionInlet > 0){
switch (BoundaryConditionInlet){
case 1:
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study);
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, Vin, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}
if (BoundaryConditionOutlet > 0){
switch (BoundaryConditionOutlet){
case 1:
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, Vout, timestep);
ScaLBL_Comm->D3Q19_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}
@ -1008,7 +1145,7 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){
}
}
void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC){
void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC, int timestep){
if (lattice_scheme.compare("D3Q7")==0){
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
@ -1029,6 +1166,30 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVe
ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
if (BoundaryConditionInlet > 0){
switch (BoundaryConditionInlet){
case 1:
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep);
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}
if (BoundaryConditionOutlet > 0){
switch (BoundaryConditionOutlet){
case 1:
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep);
break;
case 2:
Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep);
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep);
break;
}
}
ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier();
@ -1038,7 +1199,7 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVe
}
}
void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC){
void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC, int timestep){
if (lattice_scheme.compare("D3Q7")==0){
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
@ -1056,6 +1217,31 @@ void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingV
ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
if (BoundaryConditionInlet > 0){
switch (BoundaryConditionInlet){
case 1:
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
case 2:
Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep);
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
break;
}
}
if (BoundaryConditionOutlet > 0){
switch (BoundaryConditionOutlet){
case 1:
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep);
break;
case 2:
Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep);
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep);
break;
}
}
ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier();
@ -1182,6 +1368,53 @@ void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){
}
}
void ScaLBL_Poisson::WriteVis( int timestep) {
auto vis_db = db->getDatabase("Visualization");
auto format = vis_db->getWithDefault<string>( "format", "hdf5" );
DoubleArray ElectricalPotential(Nx, Ny, Nz);
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);
IO::initialize("",format,"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);
//electric potential
auto ElectricPotentialVar = std::make_shared<IO::Variable>();
//--------------------------------------------------------------------------------------------------------------------
//-------------------------------------Create Names for Variables------------------------------------------------------
if (vis_db->getWithDefault<bool>("save_electric_potential", true)) {
ElectricPotentialVar->name = "ElectricPotential";
ElectricPotentialVar->type = IO::VariableType::VolumeVariable;
ElectricPotentialVar->dim = 1;
ElectricPotentialVar->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2);
visData[0].vars.push_back(ElectricPotentialVar);
}
//--------------------------------------------------------------------------------------------------------------------
//------------------------------------Save All Variables--------------------------------------------------------------
if (vis_db->getWithDefault<bool>("save_electric_potential", true)) {
ASSERT(visData[0].vars[0]->name == "ElectricPotential");
getElectricPotential(ElectricalPotential);
Array<double> &ElectricPotentialData = visData[0].vars[0]->data;
fillData.copy(ElectricalPotential, ElectricPotentialData);
}
if (vis_db->getWithDefault<bool>("write_silo", true))
IO::writeData(timestep, visData, Dm->Comm);
//--------------------------------------------------------------------------------------------------------------------
}
//void ScaLBL_Poisson::SolveElectricField(){
// ScaLBL_Comm_Regular->SendHalo(Psi);
// ScaLBL_D3Q7_Poisson_ElectricField(NeighborList, dvcMap, dvcID, Psi, ElectricField, BoundaryConditionSolid,

View File

@ -22,6 +22,7 @@
#ifndef ScaLBL_POISSON_INC
#define ScaLBL_POISSON_INC
class ScaLBL_Poisson {
public:
ScaLBL_Poisson(int RANK, int NP, const Utilities::MPI &COMM);
@ -35,12 +36,15 @@ public:
void Create();
void Initialize(double time_conv_from_Study);
void Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study);
void Run(double *ChargeDensity, DoubleArray MembraneDistance,
bool UseSlippingVelBC, int timestep_from_Study);
void getElectricPotential(DoubleArray &ReturnValues);
void getElectricPotential_debug(int timestep);
void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y,
DoubleArray &Values_z);
void getElectricField_debug(int timestep);
void Checkpoint();
void WriteVis( int timestep);
void DummyChargeDensity(); //for debugging
@ -114,8 +118,8 @@ private:
void SolveElectricPotentialAAeven(int timestep_from_Study);
void SolveElectricPotentialAAodd(int timestep_from_Study);
//void SolveElectricField();
void SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC);
void SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC);
void SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC, int timestep);
void SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC, int timestep);
void getConvergenceLog(int timestep,double error);
double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int time_step);

View File

@ -13,8 +13,9 @@ cmake -D CMAKE_C_COMPILER:PATH=/opt/arden/openmpi/3.1.2/bin/mpicc \
-D CUDA_HOST_COMPILER="/usr/bin/gcc" \
-D USE_HDF5=1 \
-D HDF5_DIRECTORY="/opt/arden/hdf5/1.8.12" \
-D USE_DOXYGEN=true \
-D USE_CUDA=0 \
-D USE_TIMER=0 \
~/Programs/LBPM-WIA
~/Programs/LBPM
# -D HDF5_LIB="/opt/arden/hdf5/1.8.12/lib/libhdf5.a"\

View File

@ -76,11 +76,14 @@ int main(int argc, char **argv)
PoissonSolver.getElectricField_debug(timestep);
}
}
PoissonSolver.WriteVis(timestep);
}
else {
int timestep = 1;
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,false,1);
PoissonSolver.getElectricPotential_debug(1);
PoissonSolver.getElectricField_debug(1);
PoissonSolver.WriteVis(timestep);
}
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");

View File

@ -96,8 +96,10 @@ int main( int argc, char **argv )
SKIP_TIMESTEPS = 0;
if (PROTOCOL == "fractional flow")
FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
if (PROTOCOL == "seed water")
if (PROTOCOL == "seed water"){
SEED_WATER = flow_db->getWithDefault<double>( "seed_water", 0.01);
FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
}
}
/* analysis keys*/
int ANALYSIS_INTERVAL = ColorModel.timestepMax;

View File

@ -125,7 +125,7 @@ int main(int argc, char **argv)
while (timestep < Study.timestepMax){
timestep++;
PoissonSolver.Run(IonModel.ChargeDensity,SlipBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
PoissonSolver.Run(IonModel.ChargeDensity,IonModel.MembraneDistance,SlipBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
//comm.barrier();
//if (rank == 0) printf(" Poisson step %i \n",timestep);
//StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
@ -133,7 +133,7 @@ int main(int argc, char **argv)
IonModel.RunMembrane(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField,PoissonSolver.Psi); //solve for ion transport with membrane
//comm.barrier();
if (rank == 0) printf(" Membrane step %i \n",timestep);
//if (rank == 0) printf(" Membrane step %i \n",timestep);
fflush(stdout);

View File

@ -104,7 +104,6 @@ int main(int argc, char **argv)
//StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
IonModel.Run(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField); //solve for ion transport and electric potential
//if (timestep%Study.analysis_interval==0){
// Analysis.Basic(IonModel,PoissonSolver,StokesModel,timestep);
//}