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

add single cell example

Browse Source
This commit is contained in:
James McClure 2022-03-20 15:13:53 -04:00
parent 3522d35de1
commit 4927e54707
6 changed files with 180 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
common/Membrane.cpp
View File

@ -556,6 +556,9 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
membraneTag = new int [mlink];
membraneLinks = new int [2*mlink];
membraneDist = new double [2*mlink];
membraneLinkCount = mlink;
if (rank == 0) printf(" (cut %i links crossing membrane) \n",mlink);
/* construct the membrane*/
/* *
@ -744,11 +747,13 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
/* Create device copies of data structures */
ScaLBL_AllocateDeviceMemory((void **)&MembraneLinks, 2*mlink*sizeof(int));
ScaLBL_AllocateDeviceMemory((void **)&MembraneCoef, 2*mlink*sizeof(double));
ScaLBL_AllocateDeviceMemory((void **)&MembraneDistance, 2*mlink*sizeof(double));
//ScaLBL_AllocateDeviceMemory((void **)&MembraneDistance, 2*mlink*sizeof(double));
ScaLBL_AllocateDeviceMemory((void **)&MembraneDistance, Nx*Ny*Nz*sizeof(double));
ScaLBL_CopyToDevice(NeighborList, neighborList, 18*Np*sizeof(int));
ScaLBL_CopyToDevice(MembraneLinks, membraneLinks, 2*mlink*sizeof(int));
ScaLBL_CopyToDevice(MembraneDistance, membraneDist, 2*mlink*sizeof(double));
//ScaLBL_CopyToDevice(MembraneDistance, membraneDist, 2*mlink*sizeof(double));
ScaLBL_CopyToDevice(MembraneDistance, Distance.data(), Nx*Ny*Nz*sizeof(double));
if (rank == 0) printf(" Construct communication data structures... \n");
@ -1185,46 +1190,55 @@ void Membrane::RecvD3Q7AA(double *dist){
}
// std::shared_ptr<Database> db){
void Membrane::AssignCoefficients(int *Map, double *Psi, double *Distance, string method){
void Membrane::AssignCoefficients(int *Map, double *Psi, string method){
/* Assign mass transfer coefficients to the membrane data structure */
double Threshold;
double MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut;
Threshold = -55.0;
MassFractionIn = 0.0;
MassFractionOut = 0.0;
ThresholdMassFractionOut = 0.0;
ThresholdMassFractionIn = 0.0;
if (method == "Voltage Gated Potassium"){
MassFractionIn = 0.0;
MassFractionOut = 0.0;
MassFractionIn = 0.0;
ThresholdMassFractionOut = 0.0;
ThresholdMassFractionIn = -55.0;
ThresholdMassFractionIn = 1.0;
}
ScaLBL_D3Q7_Membrane_AssignLinkCoef(MembraneLinks, Map, MembraneDistance, Psi, MembraneCoef,
Threshold, MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut,
membraneLinkCount, Nx, Ny, Nz, Np);
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(-1,0,0,Map,Distance,Psi,Threshold,
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(-1,0,0,Map,MembraneDistance,Psi,Threshold,
MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut,
dvcRecvDist_X,dvcRecvLinks_X,coefficient_X,0,linkCount_X[0],recvCount_X,
Np,Nx,Ny,Nz);
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(1,0,0,Map,Distance,Psi,Threshold,
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(1,0,0,Map,MembraneDistance,Psi,Threshold,
MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut,
dvcRecvDist_x,dvcRecvLinks_x,coefficient_x,0,linkCount_x[0],recvCount_x,
Np,Nx,Ny,Nz);
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,-1,0,Map,Distance,Psi,Threshold,
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,-1,0,Map,MembraneDistance,Psi,Threshold,
MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut,
dvcRecvDist_Y,dvcRecvLinks_Y,coefficient_Y,0,linkCount_Y[0],recvCount_Y,
Np,Nx,Ny,Nz);
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,1,0,Map,Distance,Psi,Threshold,
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,1,0,Map,MembraneDistance,Psi,Threshold,
MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut,
dvcRecvDist_y,dvcRecvLinks_y,coefficient_y,0,linkCount_y[0],recvCount_y,
Np,Nx,Ny,Nz);
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,0,-1,Map,Distance,Psi,Threshold,
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,0,-1,Map,MembraneDistance,Psi,Threshold,
MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut,
dvcRecvDist_Z,dvcRecvLinks_Z,coefficient_Z,0,linkCount_Z[0],recvCount_Z,
Np,Nx,Ny,Nz);
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,0,1,Map,Distance,Psi,Threshold,
ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(0,0,1,Map,MembraneDistance,Psi,Threshold,
MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut,
dvcRecvDist_z,dvcRecvLinks_z,coefficient_z,0,linkCount_z[0],recvCount_z,
Np,Nx,Ny,Nz);

5
common/Membrane.h
View File

@ -53,6 +53,7 @@ class Membrane {
public:
int Np;
int Nx,Ny,Nz,N;
int membraneLinkCount;
int *initialNeighborList; // original neighborlist
int *NeighborList; // modified neighborlist
@ -65,7 +66,7 @@ public:
/*
* Device data structures
*/
double *MembraneLinks;
int *MembraneLinks;
double *MembraneCoef; // mass transport coefficient for the membrane
double *MembraneDistance;
@ -94,7 +95,7 @@ public:
void RecvD3Q19AA(double *dist);
void SendD3Q7AA(double *dist);
void RecvD3Q7AA(double *dist);
void AssignCoefficients(int *Map, double *Psi, double *Distance, std::string method);
void AssignCoefficients(int *Map, double *Psi, std::string method);
//......................................................................................
// Buffers to store data sent and recieved by this MPI process
double *sendbuf_x, *sendbuf_y, *sendbuf_z, *sendbuf_X, *sendbuf_Y, *sendbuf_Z;

36
example/Bubble/CreateCell.py Normal file
View File

@ -0,0 +1,36 @@
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):
C1[i,j,k] = 4.0e-6
C2[i,j,k] = 150.0e-6
dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
if (dist < 15.5 ) :
C1[i,j,k] = 140.0e-6
C2[i,j,k] = 10.0e-6
C1.tofile("cell_concentration_K_40x40x40.raw")
C2.tofile("cell_concentration_Na_40x40x40.raw")

79
example/Bubble/cell.db Normal file
View File

@ -0,0 +1,79 @@
MultiphysController {
timestepMax = 20
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_K_40x40x40.raw", "double", "cell_concentration_Na_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_InletList = 0, 0 //boundary condition for inlet; 0=periodic; 1=ion concentration; 2=ion flux
BC_OutletList = 0, 0 //boundary condition for outlet; 0=periodic; 1=ion concentration; 2=ion flux
InletValueList = 0, 0 //if ion concentration unit=[mol/m^3]; if flux (inward) unit=[mol/m^2/sec]
OutletValueList = 0, 0 //if ion concentration unit=[mol/m^3]; if flux (inward) unit=[mol/m^2/sec]
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, 1.0e-2 // 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
}

38
models/IonModel.cpp
View File

@ -629,7 +629,7 @@ void ScaLBL_IonModel::SetMembrane() {
LABEL = MembraneLabels[m];
double volume_fraction = double(label_count_global[m]) /
double((Nx - 2) * (Ny - 2) * (Nz - 2) * nprocs);
printf(" label=%d, volume fraction==%f\n", LABEL, volume_fraction);
printf(" label=%d, volume fraction = %f\n", LABEL, volume_fraction);
}
}
/* signed distance to the membrane ( - inside / + outside) */
@ -844,6 +844,7 @@ void ScaLBL_IonModel::Create() {
int neighborSize = 18 * (Np * sizeof(int));
//...........................................................................
ScaLBL_AllocateDeviceMemory((void **)&NeighborList, neighborSize);
ScaLBL_AllocateDeviceMemory((void **)&dvcMap, sizeof(int) * Np);
ScaLBL_AllocateDeviceMemory((void **)&fq,
number_ion_species * 7 * dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **)&Ci,
@ -860,6 +861,37 @@ void ScaLBL_IonModel::Create() {
if (rank == 0)
printf("LB Ion Solver: Setting up device map and neighbor list \n");
// copy the neighbor list
int *TmpMap;
TmpMap = new int[Np];
for (int k = 1; k < Nz - 1; k++) {
for (int j = 1; j < Ny - 1; j++) {
for (int i = 1; i < Nx - 1; i++) {
int idx = Map(i, j, k);
if (!(idx < 0))
TmpMap[idx] = k * Nx * Ny + j * Nx + i;
}
}
}
// check that TmpMap is valid
for (int idx = 0; idx < ScaLBL_Comm->LastExterior(); idx++) {
auto n = TmpMap[idx];
if (n > Nx * Ny * Nz) {
printf("Bad value! idx=%i \n", n);
TmpMap[idx] = Nx * Ny * Nz - 1;
}
}
for (int idx = ScaLBL_Comm->FirstInterior();
idx < ScaLBL_Comm->LastInterior(); idx++) {
auto n = TmpMap[idx];
if (n > Nx * Ny * Nz) {
printf("Bad value! idx=%i \n", n);
TmpMap[idx] = Nx * Ny * Nz - 1;
}
}
ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int) * Np);
ScaLBL_Comm->Barrier();
delete[] TmpMap;
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
comm.barrier();
@ -1263,7 +1295,9 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
//double starttime,stoptime,cputime;
//ScaLBL_Comm->Barrier(); comm.barrier();
//auto t1 = std::chrono::system_clock::now();
/* set the mass transfer coefficients for the membrane */
IonMembrane->AssignCoefficients(dvcMap, Psi, "default");
for (size_t ic = 0; ic < number_ion_species; ic++) {
timestep = 0;
while (timestep < timestepMax[ic]) {

1
models/IonModel.h
View File

@ -90,6 +90,7 @@ public:
IntArray Map;
DoubleArray Distance;
int *NeighborList;
int *dvcMap;
double *fq;
double *Ci;
double *ChargeDensity;