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

membrane starts working ok...

Browse Source
This commit is contained in:
James McClure 2022-04-08 16:44:37 -04:00
parent 25b17f996c
commit 2bb2be845a
7 changed files with 623 additions and 579 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
common/Membrane.cpp
View File

@ -420,6 +420,8 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
}
}
int Q = 7; // for D3Q7 model
/* go through the neighborlist structure */
/* count & cut the links */
if (rank == 0) printf(" Cut membrane links... \n");
@ -470,6 +472,7 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
mlink++;
}
if (Q > 7){
neighbor=Map(i-1,j-1,k);
dist=Distance(i-1,j-1,k);
if (dist*locdist < 0.0 && !(neighbor<0)){
@ -551,6 +554,7 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
}
}
}
}
/* allocate memory */
membraneTag = new int [mlink];
@ -578,8 +582,8 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
neighbor=Map(i+1,j,k);
dist=Distance(i+1,j,k);
if (dist*locdist < 0.0){
if (locdist < 0.0 && !(neighbor<0)){
if (dist*locdist < 0.0 && !(neighbor<0)){
if (locdist < 0.0 ){
localSite = 2*mlink;
neighborSite = 2*mlink+1;
}
@ -630,6 +634,8 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
mlink++;
}
if (Q > 7){
neighbor=Map(i+1,j+1,k);
dist=Distance(i+1,j+1,k);
if (dist*locdist < 0.0 && !(neighbor<0)){
@ -741,6 +747,7 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
}
}
}
}
if (rank == 0) printf(" Create device data structures... \n");
@ -1384,13 +1391,26 @@ void Membrane::AssignCoefficients(int *Map, double *Psi, string method){
ThresholdMassFractionOut = 1.0;
ThresholdMassFractionIn = 1.0;
}
if (method == "Voltage Gated Potassium"){
MassFractionIn = 0.0;
MassFractionOut = 0.0;
ThresholdMassFractionOut = 0.0;
ThresholdMassFractionIn = 1.0;
}
if (method == "impermeable"){
printf(".... impermeable membrane \n");
MassFractionIn = 0.0;
MassFractionOut = 0.0;
ThresholdMassFractionOut = 0.0;
ThresholdMassFractionIn = 0.0;
}
if (method == "Na+"){
printf(".... Na+ permeable membrane \n");
MassFractionIn = 0.5;
MassFractionOut = 0.5;
ThresholdMassFractionOut = 0.5;
ThresholdMassFractionIn = 0.5;
}
ScaLBL_D3Q7_Membrane_AssignLinkCoef(MembraneLinks, Map, MembraneDistance, Psi, MembraneCoef,
Threshold, MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut,

11
cpu/Ion.cpp
View File

@ -144,7 +144,7 @@ extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef,
iq = membrane[2*link]; ip = membrane[2*link+1];
nq = iq%Np; np = ip%Np;
fq = dist[iq]; fp = dist[ip];
fqq = (1-aq)*fq+ap*fp; fpp = (1-ap)*fp+ap*fq;
fqq = (1-aq)*fq+ap*fp; fpp = (1-ap)*fp+aq*fq;
Cq = Den[nq]; Cp = Den[np];
Cq += fqq - fq; Cp += fpp - fp;
Den[nq] = Cq; Den[np] = Cp;
@ -257,7 +257,6 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
for (n = start; n < finish; n++) {
//Load data
Ci = Den[n];
Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
@ -290,6 +289,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
f6 = dist[nr6];
// compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@ -303,6 +303,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);
@ -363,7 +365,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
for (n = start; n < finish; n++) {
//Load data
Ci = Den[n];
//Ci = Den[n];
Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
@ -383,6 +385,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
f6 = dist[5 * Np + n];
// compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@ -396,6 +399,8 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);

12
cuda/Ion.cu
View File

@ -213,7 +213,7 @@ __global__ void dvc_ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *co
iq = membrane[2*link]; ip = membrane[2*link+1];
nq = iq%Np; np = ip%Np;
fq = dist[iq]; fp = dist[ip];
fqq = (1-aq)*fq+ap*fp; fpp = (1-ap)*fp+ap*fq;
fqq = (1-aq)*fq+ap*fp; fpp = (1-ap)*fp+aq*fq;
Cq = Den[nq]; Cp = Den[np];
Cq += fqq - fq; Cp += fpp - fp;
Den[nq] = Cq; Den[np] = Cp;
@ -334,7 +334,6 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
if (n<finish) {
//Load data
Ci=Den[n];
Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
@ -367,6 +366,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
f6 = dist[nr6];
// compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@ -380,6 +380,8 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);
@ -420,6 +422,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
// q = 6
dist[nr5] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
}
}
}
@ -442,7 +445,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
if (n<finish) {
//Load data
Ci=Den[n];
//Ci = Den[n];
Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
@ -462,6 +465,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
f6 = dist[5 * Np + n];
// compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@ -475,6 +479,8 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);

37
example/SingleCell/NaCl-cell.py Normal file
View File

@ -0,0 +1,37 @@
import numpy as np
import matplotlib.pylab as plt
D=np.ones((40,40,40),dtype="uint8")
cx = 20
cy = 20
cz = 20
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 < 15.5 ) :
D[i,j,k] = 2
D.tofile("cell_40x40x40.raw")
C1=np.zeros((40,40,40),dtype="double")
C2=np.zeros((40,40,40),dtype="double")
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 < 15.5 ) :
C1[i,j,k] = 5.0e-6
C2[i,j,k] = 5.0e-6
C1.tofile("cell_concentration_Na_40x40x40.raw")
C2.tofile("cell_concentration_Cl_40x40x40.raw")

75
example/SingleCell/NaCl.db Normal file
View File

@ -0,0 +1,75 @@
MultiphysController {
timestepMax = 200
num_iter_Ion_List = 2
analysis_interval = 50
tolerance = 1.0e-9
visualization_interval = 100 // Frequency to write visualization data
analysis_interval = 50 // Frequency to perform analysis
}
Stokes {
tau = 1.0
F = 0, 0, 0
ElectricField = 0, 0, 0 //body electric field; user-input unit: [V/m]
nu_phys = 0.889e-6 //fluid kinematic viscosity; user-input unit: [m^2/sec]
}
Ions {
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
IonDiffusivityList = 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec]
IonValenceList = 1, -1 //valence charge of ions; dimensionless; positive/negative integer
IonConcentrationList = 1.0e-6, 1.0e-6 //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
}
Poisson {
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 = 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
}
Domain {
Filename = "cell_40x40x40.raw"
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz)
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
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 = true
}
Membrane {
MembraneLabels = 2
}

12
hip/Ion.hip
View File

@ -212,7 +212,7 @@ __global__ void dvc_ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *co
iq = membrane[2*link]; ip = membrane[2*link+1];
nq = iq%Np; np = ip%Np;
fq = dist[iq]; fp = dist[ip];
fqq = (1-aq)*fq+ap*fp; fpp = (1-ap)*fp+ap*fq;
fqq = (1-aq)*fq+ap*fp; fpp = (1-ap)*fp+aq*fq;
Cq = Den[nq]; Cp = Den[np];
Cq += fqq - fq; Cp += fpp - fp;
Den[nq] = Cq; Den[np] = Cp;
@ -332,7 +332,6 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
if (n<finish) {
//Load data
Ci=Den[n];
Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
@ -365,6 +364,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
f6 = dist[nr6];
// compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@ -378,6 +378,8 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);
@ -418,6 +420,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
// q = 6
dist[nr5] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
}
}
}
@ -440,7 +443,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
if (n<finish) {
//Load data
Ci=Den[n];
//Ci = Den[n];
Ex = ElectricField[n + 0 * Np];
Ey = ElectricField[n + 1 * Np];
Ez = ElectricField[n + 2 * Np];
@ -460,6 +463,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
f6 = dist[5 * Np + n];
// compute diffusive flux
Ci = f0 + f1 + f2 + f3 + f4 + f5 + f6;
flux_diffusive_x = (1.0 - 0.5 * rlx) * ((f1 - f2) - ux * Ci);
flux_diffusive_y = (1.0 - 0.5 * rlx) * ((f3 - f4) - uy * Ci);
flux_diffusive_z = (1.0 - 0.5 * rlx) * ((f5 - f6) - uz * Ci);
@ -473,6 +477,8 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
Den[n] = Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);

143
models/IonModel.cpp
View File

@ -1296,84 +1296,30 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
//ScaLBL_Comm->Barrier(); comm.barrier();
//auto t1 = std::chrono::system_clock::now();
for (size_t ic = 0; ic < number_ion_species; ic++) {
/* set the mass transfer coefficients for the membrane */
IonMembrane->AssignCoefficients(dvcMap, Psi, "default");
if (ic == 0)
IonMembrane->AssignCoefficients(dvcMap, Psi, "Na+");
else {
IonMembrane->AssignCoefficients(dvcMap, Psi, "impermeable");
}
timestep = 0;
while (timestep < timestepMax[ic]) {
//************************************************************************/
// *************ODD TIMESTEP*************//
timestep++;
//Update ion concentration and charge density
IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL
ScaLBL_D3Q7_AAodd_IonConcentration(
IonMembrane->NeighborList, &fq[ic * Np * 7], &Ci[ic * Np],
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE
/* ScaLBL_Comm->Barrier();
//--------------------------------------- Set boundary conditions -------------------------------------//
if (BoundaryConditionInlet[ic] > 0) {
switch (BoundaryConditionInlet[ic]) {
case 1:
ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], timestep);
break;
case 21:
ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
&Velocity[2 * Np], timestep);
break;
case 22:
ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
&Velocity[2 * Np], timestep);
break;
case 23:
ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
&Velocity[2 * Np], &ElectricField[2 * Np],
IonDiffusivity[ic], IonValence[ic], Vt, timestep);
break;
}
}
if (BoundaryConditionOutlet[ic] > 0) {
switch (BoundaryConditionOutlet[ic]) {
case 1:
ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], timestep);
break;
case 21:
ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
&Velocity[2 * Np], timestep);
break;
case 22:
ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
&Velocity[2 * Np], timestep);
break;
case 23:
ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
&Velocity[2 * Np], &ElectricField[2 * Np],
IonDiffusivity[ic], IonValence[ic], Vt, timestep);
break;
}
}
*/
//----------------------------------------------------------------------------------------------------//
ScaLBL_D3Q7_AAodd_IonConcentration(IonMembrane->NeighborList, &fq[ic * Np * 7],
&Ci[ic * Np], 0,
ScaLBL_Comm->LastExterior(), Np);
//LB-Ion collison
IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL
ScaLBL_D3Q7_AAodd_Ion(
IonMembrane->NeighborList, &fq[ic * Np * 7], &Ci[ic * Np],
&FluxDiffusive[3 * ic * Np], &FluxAdvective[3 * ic * Np],
&FluxElectrical[3 * ic * Np], Velocity, ElectricField,
IonDiffusivity[ic], IonValence[ic], rlx[ic], Vt,
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAodd_Ion(
IonMembrane->NeighborList, &fq[ic * Np * 7], &Ci[ic * Np],
&FluxDiffusive[3 * ic * Np], &FluxAdvective[3 * ic * Np],
@ -1381,6 +1327,9 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
IonDiffusivity[ic], IonValence[ic], rlx[ic], Vt, 0,
ScaLBL_Comm->LastExterior(), Np);
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);
@ -1390,81 +1339,27 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
// *************EVEN TIMESTEP*************//
timestep++;
//Update ion concentration and charge density
IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL
ScaLBL_D3Q7_AAeven_IonConcentration(
&fq[ic * Np * 7], &Ci[ic * Np], ScaLBL_Comm->FirstInterior(),
ScaLBL_Comm->LastInterior(), Np);
IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE
ScaLBL_Comm->Barrier();
//--------------------------------------- Set boundary conditions -------------------------------------//
/*if (BoundaryConditionInlet[ic] > 0) {
switch (BoundaryConditionInlet[ic]) {
case 1:
ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], timestep);
break;
case 21:
ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],q
&Velocity[2 * Np], timestep);
break;
case 22:
ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
&Velocity[2 * Np], timestep);
break;
case 23:
ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cin[ic], tau[ic],
&Velocity[2 * Np], &ElectricField[2 * Np],
IonDiffusivity[ic], IonValence[ic], Vt, timestep);
break;
}
}
if (BoundaryConditionOutlet[ic] > 0) {
switch (BoundaryConditionOutlet[ic]) {
case 1:
ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], timestep);
break;
case 21:
ScaLBL_Comm->D3Q7_Ion_Flux_Diff_BC_Z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
&Velocity[2 * Np], timestep);
break;
case 22:
ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvc_BC_Z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
&Velocity[2 * Np], timestep);
break;
case 23:
ScaLBL_Comm->D3Q7_Ion_Flux_DiffAdvcElec_BC_Z(
IonMembrane->NeighborList, &fq[ic * Np * 7], Cout[ic], tau[ic],
&Velocity[2 * Np], &ElectricField[2 * Np],
IonDiffusivity[ic], IonValence[ic], Vt, timestep);
break;
}
}
*/
//----------------------------------------------------------------------------------------------------//
ScaLBL_D3Q7_AAeven_IonConcentration(&fq[ic * Np * 7], &Ci[ic * Np],
0, ScaLBL_Comm->LastExterior(),
Np);
//LB-Ion collison
IonMembrane->SendD3Q7AA(&fq[ic * Np * 7]); //READ FORM NORMAL
ScaLBL_D3Q7_AAeven_Ion(
&fq[ic * Np * 7], &Ci[ic * Np], &FluxDiffusive[3 * ic * Np],
&FluxAdvective[3 * ic * Np], &FluxElectrical[3 * ic * Np],
Velocity, ElectricField, IonDiffusivity[ic], IonValence[ic],
rlx[ic], Vt, ScaLBL_Comm->FirstInterior(),
ScaLBL_Comm->LastInterior(), Np);
IonMembrane->RecvD3Q7AA(&fq[ic * Np * 7]); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAeven_Ion(
&fq[ic * Np * 7], &Ci[ic * Np], &FluxDiffusive[3 * ic * Np],
&FluxAdvective[3 * ic * Np], &FluxElectrical[3 * ic * Np],
Velocity, ElectricField, IonDiffusivity[ic], IonValence[ic],
rlx[ic], Vt, 0, ScaLBL_Comm->LastExterior(), Np);
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);