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

Merge branch 'membrane' into tmp

Browse Source
This commit is contained in:
James E McClure 2022-03-27 16:03:26 -04:00
commit f99af45c9b
19 changed files with 3332 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
common/Domain.cpp
View File

@ -1543,7 +1543,7 @@ void Domain::ReadFromFile(const std::string &Filename,
} else {
// Recieve the subdomain from rank = 0
//printf("Ready to recieve data %i at process %i \n", N,rank);
Comm.recv(id.data(), N, 0, 15);
Comm.recv(UserData, N, 0, 15);
}
Comm.barrier();
}

1415
common/Membrane.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

183
common/Membrane.h Normal file
View File

@ -0,0 +1,183 @@
/* Flow adaptor class for multiphase flow methods */
#ifndef ScaLBL_Membrane_INC
#define ScaLBL_Membrane_INC
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include "common/ScaLBL.h"
/**
* \brief Unpack D3Q19 distributions after communication using links determined based on membrane location
* @param q - index for distribution based on D3Q19 discrete velocity structure
* @param list - list of distributions to communicate
* @param links - list of active links based on the membrane location
* @param start - index to start parsing the list
* @param count - number of values to unppack
* @param recvbuf - memory buffer where recieved values have been stored
* @param dist - memory buffer to hold the distributions
* @param N - size of the distributions (derived from Domain structure)
*/
extern "C" void Membrane_D3Q19_Unpack(int q, int *list, int *links, int start, int count, double *recvbuf, double *dist, int N);
/**
* \brief Set custom link rules for D3Q19 distribution based on membrane location
* @param q - index for distribution based on D3Q19 discrete velocity structure
* @param list - list of distributions to communicate
* @param links - list of active links based on the membrane location
* @param coef - coefficient to determine the local mass transport for each membrane link
* @param start - index to start parsing the list
* @param offset - offset to start reading membrane links
* @param count - number of values to unppack
* @param recvbuf - memory buffer where recieved values have been stored
* @param dist - memory buffer to hold the distributions
* @param N - size of the distributions (derived from Domain structure)
*/
extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links, double *coef, int start, int offset,
int linkCount, double *recvbuf, double *dist, int N);
/**
* \class Membrane
* @brief
* The Membrane class operates on ScaLBL data structures to insert membrane
*
*/
class Membrane {
public:
int Np;
int Nx,Ny,Nz,N;
int membraneLinkCount;
int *initialNeighborList; // original neighborlist
int *NeighborList; // modified neighborlist
/* host data structures */
int *membraneLinks; // D3Q19 links that cross membrane
int *membraneTag; // label each link in the membrane
double *membraneDist; // distance to membrane for each linked site
/*
* Device data structures
*/
int *MembraneLinks;
double *MembraneCoef; // mass transport coefficient for the membrane
double *MembraneDistance;
/**
* \brief Create a flow adaptor to operate on the LB model
* @param ScaLBL - originating data structures
* @param neighborList - list of neighbors for each site
*/
Membrane(std::shared_ptr <Domain> Dm, int *initialNeighborList, int Nsites);
/**
* \brief Destructor
*/
~Membrane();
/**
* \brief Create membrane
* \details Create membrane structure from signed distance function
* @param Dm - domain structure
* @param Distance - signed distance to membrane
* @param Map - mapping between regular layout and compact layout
*/
int Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArray &Map);
void SendD3Q19AA(double *dist);
void RecvD3Q19AA(double *dist);
void SendD3Q7AA(double *dist);
void RecvD3Q7AA(double *dist);
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;
double *sendbuf_xy, *sendbuf_yz, *sendbuf_xz, *sendbuf_Xy, *sendbuf_Yz, *sendbuf_xZ;
double *sendbuf_xY, *sendbuf_yZ, *sendbuf_Xz, *sendbuf_XY, *sendbuf_YZ, *sendbuf_XZ;
double *recvbuf_x, *recvbuf_y, *recvbuf_z, *recvbuf_X, *recvbuf_Y, *recvbuf_Z;
double *recvbuf_xy, *recvbuf_yz, *recvbuf_xz, *recvbuf_Xy, *recvbuf_Yz, *recvbuf_xZ;
double *recvbuf_xY, *recvbuf_yZ, *recvbuf_Xz, *recvbuf_XY, *recvbuf_YZ, *recvbuf_XZ;
//......................................................................................
private:
bool Lock; // use Lock to make sure only one call at a time to protect data in transit
int sendtag, recvtag;
int iproc,jproc,kproc;
int nprocx,nprocy,nprocz;
// Give the object it's own MPI communicator
RankInfoStruct rank_info;
Utilities::MPI MPI_COMM_SCALBL; // MPI Communicator for this domain
MPI_Request req1[18],req2[18];
/**
* \brief Set up membrane communication
* \details associate p2p communication links to membrane where necessary
* returns the number of membrane links
* regular communications are stored in the first part of the list
* membrane communications are stored in the last part of the list
* @param Cqx - discrete velocity (x)
* @param Cqy - discrete velocity (y)
* @param Cqz - discrete velocity (z)
* @param list - list of recieved values
* @param count - number recieved values
* @param d3q19_recvlist - device array with the saved list
* @param d3q19_linkList - sorted list with regular and membrane links
* @param Distance - signed distance to membrane
* @param Map - data structure used to define mapping between dense and sparse representation
* */
int D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, const int *list, int start, int count,
int *d3q19_recvlist, int *d3q19_linkList, DoubleArray &Distance, IntArray &Map);
//......................................................................................
// MPI ranks for all 18 neighbors
//......................................................................................
// These variables are all private to prevent external things from modifying them!!
//......................................................................................
int rank;
int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z;
int rank_xy,rank_XY,rank_xY,rank_Xy;
int rank_xz,rank_XZ,rank_xZ,rank_Xz;
int rank_yz,rank_YZ,rank_yZ,rank_Yz;
//......................................................................................
int SendCount, RecvCount, CommunicationCount;
//......................................................................................
int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z;
int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ;
int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ;
//......................................................................................
int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z;
int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ;
int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ;
//......................................................................................
int linkCount_x[5], linkCount_y[5], linkCount_z[5], linkCount_X[5], linkCount_Y[5], linkCount_Z[5];
int linkCount_xy, linkCount_yz, linkCount_xz, linkCount_Xy, linkCount_Yz, linkCount_xZ;
int linkCount_xY, linkCount_yZ, linkCount_Xz, linkCount_XY, linkCount_YZ, linkCount_XZ;
//......................................................................................
// Send buffers that reside on the compute device
int *dvcSendList_x, *dvcSendList_y, *dvcSendList_z, *dvcSendList_X, *dvcSendList_Y, *dvcSendList_Z;
int *dvcSendList_xy, *dvcSendList_yz, *dvcSendList_xz, *dvcSendList_Xy, *dvcSendList_Yz, *dvcSendList_xZ;
int *dvcSendList_xY, *dvcSendList_yZ, *dvcSendList_Xz, *dvcSendList_XY, *dvcSendList_YZ, *dvcSendList_XZ;
// Recieve buffers that reside on the compute device
int *dvcRecvList_x, *dvcRecvList_y, *dvcRecvList_z, *dvcRecvList_X, *dvcRecvList_Y, *dvcRecvList_Z;
int *dvcRecvList_xy, *dvcRecvList_yz, *dvcRecvList_xz, *dvcRecvList_Xy, *dvcRecvList_Yz, *dvcRecvList_xZ;
int *dvcRecvList_xY, *dvcRecvList_yZ, *dvcRecvList_Xz, *dvcRecvList_XY, *dvcRecvList_YZ, *dvcRecvList_XZ;
// Link lists that reside on the compute device
int *dvcRecvLinks_x, *dvcRecvLinks_y, *dvcRecvLinks_z, *dvcRecvLinks_X, *dvcRecvLinks_Y, *dvcRecvLinks_Z;
int *dvcRecvLinks_xy, *dvcRecvLinks_yz, *dvcRecvLinks_xz, *dvcRecvLinks_Xy, *dvcRecvLinks_Yz, *dvcRecvLinks_xZ;
int *dvcRecvLinks_xY, *dvcRecvLinks_yZ, *dvcRecvLinks_Xz, *dvcRecvLinks_XY, *dvcRecvLinks_YZ, *dvcRecvLinks_XZ;
// Recieve buffers for the distributions
int *dvcRecvDist_x, *dvcRecvDist_y, *dvcRecvDist_z, *dvcRecvDist_X, *dvcRecvDist_Y, *dvcRecvDist_Z;
int *dvcRecvDist_xy, *dvcRecvDist_yz, *dvcRecvDist_xz, *dvcRecvDist_Xy, *dvcRecvDist_Yz, *dvcRecvDist_xZ;
int *dvcRecvDist_xY, *dvcRecvDist_yZ, *dvcRecvDist_Xz, *dvcRecvDist_XY, *dvcRecvDist_YZ, *dvcRecvDist_XZ;
//......................................................................................
// mass transfer coefficient arrays
double *coefficient_x, *coefficient_X, *coefficient_y, *coefficient_Y, *coefficient_z, *coefficient_Z;
//......................................................................................
};
#endif

9
common/ScaLBL.cpp
View File

@ -1,5 +1,5 @@
/*
Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University
Copyright 2013--2022 James E. McClure, Virginia Polytechnic & State University
Copyright Equnior ASA
This file is part of the Open Porous Media project (OPM).
@ -15,10 +15,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "common/ScaLBL.h"
#include <chrono>
ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr <Domain> Dm){
//......................................................................................
Lock=false; // unlock the communicator
@ -366,7 +364,6 @@ ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr <Domain> Dm){
}
ScaLBL_Communicator::~ScaLBL_Communicator()
{
@ -888,8 +885,6 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis
idx=Map(n);
//if (rank == 0) printf("r: mapped n=%d\n",idx);
TempBuffer[i]=idx;
}
ScaLBL_CopyToDevice(dvcRecvDist_x,TempBuffer,5*recvCount_x*sizeof(int));
@ -1046,7 +1041,6 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis
return(Np);
}
void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC)
{
@ -1448,7 +1442,6 @@ void ScaLBL_Communicator::SolidSlippingVelocityBCD3Q19(double *fq, double *zeta_
void ScaLBL_Communicator::SendD3Q19AA(double *dist){
// NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2
if (Lock==true){
ERROR("ScaLBL Error (SendD3Q19): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?");
}

28
common/ScaLBL.h
View File

@ -217,6 +217,25 @@ extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int
*/
extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz);
// MEMBRANE MODEL
extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef, double *dist, double *Den, int memLinks, int Np);
extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef,
double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int memLinks, int Nx, int Ny, int Nz, int Np);
extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(
const int Cqx, const int Cqy, int const Cqz,
int *Map, double *Distance, double *Psi, double Threshold,
double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count,
const int N, const int Nx, const int Ny, const int Nz);
extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q,
int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count,
double *recvbuf, double *dist, int N, double *coef);
// GREYSCALE MODEL (Single-component)
extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *Dist, int Np, double Den);
@ -702,6 +721,14 @@ public:
double GetPerformance(int *NeighborList, double *fq, int Np);
int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np, int width);
/**
* \brief Create membrane data structure
* - cut lattice links based on distance map
* @param Distance - signed distance to membrane
* @param neighborList - data structure that retains lattice links
* @param Np - number of lattice sites
* @param width - halo width for the model
*/
void Barrier(){
ScaLBL_DeviceBarrier();
MPI_COMM_SCALBL.barrier();
@ -782,7 +809,6 @@ private:
int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ;
int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ;
//......................................................................................
int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z;
int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ;
int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ;

1
cpu/D3Q19.cpp
View File

@ -48,6 +48,7 @@ extern "C" void ScaLBL_D3Q19_Unpack(int q, int *list, int start, int count,
}
}
extern "C" void ScaLBL_D3Q19_AA_Init(double *f_even, double *f_odd, int Np) {
int n;
for (n = 0; n < Np; n++) {

207
cpu/Ion.cpp
View File

@ -1,4 +1,156 @@
#include <stdio.h>
#include <math.h>
extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef,
double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int memLinks, int Nx, int Ny, int Nz, int Np){
int link,iq,ip,nq,np,nqm,npm;
double aq, ap, membranePotential;
/* Interior Links */
for (link=0; link<memLinks; link++){
// inside //outside
aq = MassFractionIn; ap = MassFractionOut;
iq = membrane[2*link]; ip = membrane[2*link+1];
nq = iq%Np; np = ip%Np;
nqm = Map[nq]; npm = Map[np]; // strided layout
/* membrane potential for this link */
membranePotential = Psi[nqm] - Psi[npm];
if (membranePotential > Threshold){
aq = ThresholdMassFractionIn; ap = ThresholdMassFractionOut;
}
/* Save the mass transfer coefficients */
coef[2*link] = aq; coef[2*link+1] = ap;
}
}
extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(
const int Cqx, const int Cqy, int const Cqz,
int *Map, double *Distance, double *Psi, double Threshold,
double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count,
const int N, const int Nx, const int Ny, const int Nz) {
//....................................................................................
// Unack distribution from the recv buffer
// Distribution q matche Cqx, Cqy, Cqz
// swap rule means that the distributions in recvbuf are OPPOSITE of q
// dist may be even or odd distributions stored by stream layout
//....................................................................................
int n, idx, link, nqm, npm, i, j, k;
double distanceLocal, distanceNonlocal;
double psiLocal, psiNonlocal, membranePotential;
double ap,aq; // coefficient
/* second enforce custom rule for membrane links */
for (link = nlinks; link < count; link++) {
// get the index for the recv list (deal with reordering of links)
idx = d3q7_linkList[link]; // THINK start NEEDS TO BE HERE
// get the distribution index
n = d3q7_recvlist[start+idx];
// get the index in strided layout
nqm = Map[n];
distanceLocal = Distance[nqm];
psiLocal = Psi[nqm];
// Get the 3-D indices from the send process
k = nqm/(Nx*Ny); j = (nqm-Nx*Ny*k)/Nx; i = nqm-Nx*Ny*k-Nx*j;
// Streaming link the non-local distribution
i -= Cqx; j -= Cqy; k -= Cqz;
npm = k*Nx*Ny + j*Nx + i;
distanceNonlocal = Distance[npm];
psiNonlocal = Psi[npm];
membranePotential = psiLocal - psiNonlocal;
aq = MassFractionIn;
ap = MassFractionOut;
/* link is inside membrane */
if (distanceLocal > 0.0){
if (membranePotential < Threshold*(-1.0)){
ap = MassFractionIn;
aq = MassFractionOut;
}
else {
ap = ThresholdMassFractionIn;
aq = ThresholdMassFractionOut;
}
}
else if (membranePotential > Threshold){
aq = ThresholdMassFractionIn;
ap = ThresholdMassFractionOut;
}
// update link based on mass transfer coefficients
coef[2*(link-nlinks)] = aq;
coef[2*(link-nlinks)+1] = ap;
}
}
extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q,
int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count,
double *recvbuf, double *dist, int N, double *coef) {
//....................................................................................
// Unack distribution from the recv buffer
// Distribution q matche Cqx, Cqy, Cqz
// swap rule means that the distributions in recvbuf are OPPOSITE of q
// dist may be even or odd distributions stored by stream layout
//....................................................................................
int n, idx, link;
double fq,fp,fqq,ap,aq; // coefficient
/* First unpack the regular links */
for (link = 0; link < nlinks; link++) {
// get the index for the recv list (deal with reordering of links)
idx = d3q7_linkList[link];
// get the distribution index
n = d3q7_recvlist[start+idx];
if (!(n < 0)){
fp = recvbuf[start + idx];
dist[q * N + n] = fp;
}
//printf(" site=%i, index=%i, value = %e \n",n,idx,fp);
}
/* second enforce custom rule for membrane links */
for (link = nlinks; link < count; link++) {
// get the index for the recv list (deal with reordering of links)
idx = d3q7_linkList[link];
// get the distribution index
n = d3q7_recvlist[start+idx];
// update link based on mass transfer coefficients
if (!(n < 0)){
aq = coef[2*(link-nlinks)];
ap = coef[2*(link-nlinks)+1];
fq = dist[q * N + n];
fp = recvbuf[start + idx];
fqq = (1-aq)*fq+ap*fp;
dist[q * N + n] = fqq;
}
//printf(" LINK: site=%i, index=%i \n", n, idx);
}
}
extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef,
double *dist, double *Den, int memLinks, int Np){
int link,iq,ip,nq,np;
double aq, ap, fq, fp, fqq, fpp, Cq, Cp;
for (link=0; link<memLinks; link++){
// inside //outside
aq = coef[2*link]; ap = coef[2*link+1];
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;
Cq = Den[nq]; Cp = Den[np];
Cq += fqq - fq; Cp += fpp - fp;
Den[nq] = Cq; Den[np] = Cp;
dist[iq] = fqq; dist[ip] = fpp;
}
}
extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList,
double *dist, double *Den,
@ -99,6 +251,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
double Ex, Ey, Ez; //electrical field
double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z;
double f0, f1, f2, f3, f4, f5, f6;
double X,Y,Z,factor_x, factor_y, factor_z;
int nr1, nr2, nr3, nr4, nr5, nr6;
for (n = start; n < finish; n++) {
@ -149,33 +302,48 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist,
FluxElectrical[n + 0 * Np] = uEPx * Ci;
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);
Z = 4.0 * (uz + uEPz);
factor_x = X / sqrt(1 + X*X);
factor_y = Y / sqrt(1 + Y*Y);
factor_z = Z / sqrt(1 + Z*Z);
// q=0
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1
dist[nr2] =
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
//f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
// q=2
dist[nr1] =
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
//f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
// q = 3
dist[nr4] =
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y );
//f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
// q = 4
dist[nr3] =
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
//f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
// q = 5
dist[nr6] =
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
//f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
// q = 6
dist[nr5] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
}
}
@ -190,6 +358,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
double Ex, Ey, Ez; //electrical field
double flux_diffusive_x, flux_diffusive_y, flux_diffusive_z;
double f0, f1, f2, f3, f4, f5, f6;
double X,Y,Z, factor_x, factor_y, factor_z;
for (n = start; n < finish; n++) {
@ -226,33 +395,47 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(
FluxElectrical[n + 0 * Np] = uEPx * Ci;
FluxElectrical[n + 1 * Np] = uEPy * Ci;
FluxElectrical[n + 2 * Np] = uEPz * Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);
Z = 4.0 * (uz + uEPz);
factor_x = X / sqrt(1 + X*X);
factor_y = Y / sqrt(1 + Y*Y);
factor_z = Z / sqrt(1 + Z*Z);
// q=0
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1
dist[1 * Np + n] =
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
//f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
// q=2
dist[2 * Np + n] =
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
//f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
// q = 3
dist[3 * Np + n] =
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y);
//f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
// q = 4
dist[4 * Np + n] =
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
//f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
// q = 5
dist[5 * Np + n] =
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
//f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
// q = 6
dist[6 * Np + n] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
//f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
}
}

42
cpu/MembraneHelper.cpp Normal file
View File

@ -0,0 +1,42 @@
extern "C" void Membrane_D3Q19_Unpack(int q, int *list, int *links, int start, int linkCount,
double *recvbuf, double *dist, int N) {
//....................................................................................
// Unack distribution from the recv buffer
// Distribution q matche Cqx, Cqy, Cqz
// swap rule means that the distributions in recvbuf are OPPOSITE of q
// dist may be even or odd distributions stored by stream layout
//....................................................................................
int n, idx, link;
for (link=0; link<linkCount; link++){
idx = links[start+link];
// Get the value from the list -- note that n is the index is from the send (non-local) process
n = list[start + idx];
// unpack the distribution to the proper location
if (!(n < 0))
dist[q * N + n] = recvbuf[start + idx];
}
}
extern "C" void Membrane_D3Q19_Transport(int q, int *list, int *links, double *coef, int start, int offset,
int linkCount, double *recvbuf, double *dist, int N){
//....................................................................................
// Unack distribution from the recv buffer
// Distribution q matche Cqx, Cqy, Cqz
// swap rule means that the distributions in recvbuf are OPPOSITE of q
// dist may be even or odd distributions stored by stream layout
//....................................................................................
int n, idx, link;
double alpha;
for (link=offset; link<linkCount; link++){
idx = list[start+link];
// Get the value from the list -- note that n is the index is from the send (non-local) process
n = list[start + idx];
alpha = coef[start + idx];
// unpack the distribution to the proper location
if (!(n < 0))
dist[q * N + n] = alpha*recvbuf[start + idx];
}
}

318
cuda/Ion.cu
View File

@ -5,6 +5,182 @@
#define NBLOCKS 1024
#define NTHREADS 256
__global__ void dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef,
double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int memLinks, int Nx, int Ny, int Nz, int Np){
int link,iq,ip,nq,np,nqm,npm;
double aq, ap, membranePotential;
/* Interior Links */
int S = memLinks/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
link = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
if (link < memLinks) {
// inside //outside
aq = MassFractionIn; ap = MassFractionOut;
iq = membrane[2*link]; ip = membrane[2*link+1];
nq = iq%Np; np = ip%Np;
nqm = Map[nq]; npm = Map[np]; // strided layout
/* membrane potential for this link */
membranePotential = Psi[nqm] - Psi[npm];
if (membranePotential > Threshold){
aq = ThresholdMassFractionIn; ap = ThresholdMassFractionOut;
}
/* Save the mass transfer coefficients */
coef[2*link] = aq; coef[2*link+1] = ap;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(
const int Cqx, const int Cqy, int const Cqz,
int *Map, double *Distance, double *Psi, double Threshold,
double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count,
const int N, const int Nx, const int Ny, const int Nz) {
//....................................................................................
// Unack distribution from the recv buffer
// Distribution q matche Cqx, Cqy, Cqz
// swap rule means that the distributions in recvbuf are OPPOSITE of q
// dist may be even or odd distributions stored by stream layout
//....................................................................................
int n, idx, link, nqm, npm, i, j, k;
double distanceLocal, distanceNonlocal;
double psiLocal, psiNonlocal, membranePotential;
double ap,aq; // coefficient
/* second enforce custom rule for membrane links */
int S = (count-nlinks)/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
link = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + nlinks;
if (link < count) {
// get the index for the recv list (deal with reordering of links)
idx = d3q7_linkList[link]; // THINK start NEEDS TO BE HERE
// get the distribution index
n = d3q7_recvlist[start+idx];
// get the index in strided layout
nqm = Map[n];
distanceLocal = Distance[nqm];
psiLocal = Psi[nqm];
// Get the 3-D indices from the send process
k = nqm/(Nx*Ny); j = (nqm-Nx*Ny*k)/Nx; i = nqm-Nx*Ny*k-Nx*j;
// Streaming link the non-local distribution
i -= Cqx; j -= Cqy; k -= Cqz;
npm = k*Nx*Ny + j*Nx + i;
distanceNonlocal = Distance[npm];
psiNonlocal = Psi[npm];
membranePotential = psiLocal - psiNonlocal;
aq = MassFractionIn;
ap = MassFractionOut;
/* link is inside membrane */
if (distanceLocal > 0.0){
if (membranePotential < Threshold*(-1.0)){
ap = MassFractionIn;
aq = MassFractionOut;
}
else {
ap = ThresholdMassFractionIn;
aq = ThresholdMassFractionOut;
}
}
else if (membranePotential > Threshold){
aq = ThresholdMassFractionIn;
ap = ThresholdMassFractionOut;
}
// update link based on mass transfer coefficients
coef[2*(link-nlinks)] = aq;
coef[2*(link-nlinks)+1] = ap;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_Membrane_Unpack(int q,
int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count,
double *recvbuf, double *dist, int N, double *coef) {
//....................................................................................
// Unack distribution from the recv buffer
// Distribution q matche Cqx, Cqy, Cqz
// swap rule means that the distributions in recvbuf are OPPOSITE of q
// dist may be even or odd distributions stored by stream layout
//....................................................................................
int n, idx, link;
double fq,fp,fqq,ap,aq; // coefficient
/* second enforce custom rule for membrane links */
int S = count/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
link = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
/* First unpack the regular links */
if (link < nlinks) {
// get the index for the recv list (deal with reordering of links)
idx = d3q7_linkList[link];
// get the distribution index
n = d3q7_recvlist[start+idx];
if (!(n < 0)){
fp = recvbuf[start + idx];
dist[q * N + n] = fp;
}
}
else if (link < count){
/* second enforce custom rule for membrane links */
// get the index for the recv list (deal with reordering of links)
idx = d3q7_linkList[link];
// get the distribution index
n = d3q7_recvlist[start+idx];
// update link based on mass transfer coefficients
if (!(n < 0)){
aq = coef[2*(link-nlinks)];
ap = coef[2*(link-nlinks)+1];
fq = dist[q * N + n];
fp = recvbuf[start + idx];
fqq = (1-aq)*fq+ap*fp;
dist[q * N + n] = fqq;
}
}
}
}
__global__ void dvc_ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef,
double *dist, double *Den, int memLinks, int Np){
int link,iq,ip,nq,np;
double aq, ap, fq, fp, fqq, fpp, Cq, Cp;
int S = memLinks/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
link = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
if (link < memLinks){
// inside //outside
aq = coef[2*link]; ap = coef[2*link+1];
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;
Cq = Den[nq]; Cp = Den[np];
Cq += fqq - fq; Cp += fpp - fp;
Den[nq] = Cq; Den[np] = Cp;
dist[iq] = fqq; dist[ip] = fpp;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *dist, double *Den, int start, int finish, int Np){
int n,nread;
double fq,Ci;
@ -106,6 +282,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
double X,Y,Z,factor_x,factor_y,factor_z;
int nr1,nr2,nr3,nr4,nr5,nr6;
int S = Np/NBLOCKS/NTHREADS + 1;
@ -160,34 +337,47 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);
Z = 4.0 * (uz + uEPz);
factor_x = X / sqrt(1 + X*X);
factor_y = Y / sqrt(1 + Y*Y);
factor_z = Z / sqrt(1 + Z*Z);
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
//dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1
dist[nr2] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx));
//dist[nr2] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr2] =
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
//f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
// q=2
dist[nr1] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx));
//dist[nr1] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr1] =
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
//f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
// q = 3
dist[nr4] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy));
//dist[nr4] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr4] =
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y );
//f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
// q = 4
dist[nr3] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy));
//dist[nr3] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr3] =
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
//f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
// q = 5
dist[nr6] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz));
//dist[nr6] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr6] =
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
//f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
// q = 6
dist[nr5] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz));
//dist[nr5] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr5] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
}
}
}
@ -201,6 +391,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
double X,Y,Z,factor_x,factor_y,factor_z;
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
@ -242,33 +433,46 @@ __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;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);
Z = 4.0 * (uz + uEPz);
factor_x = X / sqrt(1 + X*X);
factor_y = Y / sqrt(1 + Y*Y);
factor_z = Z / sqrt(1 + Z*Z);
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
//dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1
dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx));
//dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[1 * Np + n] =
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
//f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
// q=2
dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx));
//dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[2 * Np + n] =
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
//f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
// q = 3
dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy));
//dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[3 * Np + n] =
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y);
//f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
// q = 4
dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy));
//dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[4 * Np + n] =
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
//f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
// q = 5
dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz));
//dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[5 * Np + n] =
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
//f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
// q = 6
dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz));
//dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[6 * Np + n] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
//f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
}
}
}
@ -419,3 +623,61 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity
}
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef,
double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int memLinks, int Nx, int Ny, int Nz, int Np){
dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef<<<NBLOCKS,NTHREADS >>>(membrane, Map, Distance, Psi, coef,
Threshold, MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut,
memLinks, Nx, Ny, Nz, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef: %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(
const int Cqx, const int Cqy, int const Cqz,
int *Map, double *Distance, double *Psi, double Threshold,
double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count,
const int N, const int Nx, const int Ny, const int Nz) {
dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo<<<NBLOCKS,NTHREADS >>>(
Cqx, Cqy, Cqz, Map, Distance, Psi, Threshold,
MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut,
d3q7_recvlist, d3q7_linkList, coef, start, nlinks, count, N, Nx, Ny, Nz);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo: %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q,
int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count,
double *recvbuf, double *dist, int N, double *coef) {
dvc_ScaLBL_D3Q7_Membrane_Unpack<<<NBLOCKS,NTHREADS >>>(q, d3q7_recvlist, d3q7_linkList, start, nlinks, count,
recvbuf, dist, N, coef) ;
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_Unpack: %s \n",cudaGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef,
double *dist, double *Den, int memLinks, int Np){
dvc_ScaLBL_D3Q7_Membrane_IonTransport<<<NBLOCKS,NTHREADS >>>(membrane, coef, dist, Den, memLinks, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_IonTransport: %s \n",cudaGetErrorString(err));
}
}

18
example/Bubble/CreateBubble.py Normal file
View File

@ -0,0 +1,18 @@
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 < 12.5 ) :
D[i,j,k] = 2
D.tofile("bubble_40x40x40.raw")

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

@ -0,0 +1,77 @@
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")
C3=np.zeros((40,40,40),dtype="double")
C4=np.zeros((40,40,40),dtype="double")
C5=np.zeros((40,40,40),dtype="double")
C6=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] = 4.0e-6 # K
C2[i,j,k] = 150.0e-6 # Na
C3[i,j,k] = 116.0e-6 # Cl
C4[i,j,k] = 29.0e-6 # HC03
#C5[i,j,k] = 2.4e-6 # Ca
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] = 145.0e-6
C2[i,j,k] = 12.0e-6
C3[i,j,k] = 4.0e-6
C4[i,j,k] = 12.0e-6 # 12 mmol / L
#C5[i,j,k] = 0.10e-6 # 100 nmol / L
# add up the total charge to make sure it is zero
TotalCharge = 0
for i in range(0, 40):
for j in range (0, 40):
for k in range (0,40):
TotalCharge += C1[i,j,k] + C2[i,j,k] - C3[i,j,k] - C4[i,j,k]
TotalCharge /= (40*40*40)
print("Total charge " + str(TotalCharge))
for i in range(0, 40):
for j in range (0, 40):
for k in range (0,40):
if TotalCharge < 0 :
# need more cation
C5[i,j,k] = abs(TotalCharge)
C6[i,j,k] = 0.0
else :
# need more anion
C5[i,j,k] = 0.0
C6[i,j,k] = abs(TotalCharge)
C1.tofile("cell_concentration_K_40x40x40.raw")
C2.tofile("cell_concentration_Na_40x40x40.raw")
C3.tofile("cell_concentration_Cl_40x40x40.raw")
C4.tofile("cell_concentration_HCO3_40x40x40.raw")
C5.tofile("cell_concentration_cation_40x40x40.raw")
C6.tofile("cell_concentration_anion_40x40x40.raw")

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

@ -0,0 +1,75 @@
MultiphysController {
timestepMax = 60
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", "cell_concentration_Cl_40x40x40.raw", "double", "cell_concentration_HCO3_40x40x40.raw", "double", "cell_concentration_anion_40x40x40.raw", "double", "cell_concentration_cation_40x40x40.raw", "double"
temperature = 293.15 //unit [K]
number_ion_species = 6 //number of ions
tauList = 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
IonDiffusivityList = 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9 //user-input unit: [m^2/sec]
IonValenceList = 1, 1, -1, -1, 1, -1 //valence charge of ions; dimensionless; positive/negative integer
IonConcentrationList = 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 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, 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
}

314
hip/Ion.hip
View File

@ -5,6 +5,179 @@
#define NBLOCKS 1024
#define NTHREADS 256
__global__ void dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef,
double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int memLinks, int Nx, int Ny, int Nz, int Np){
int link,iq,ip,nq,np,nqm,npm;
double aq, ap, membranePotential;
/* Interior Links */
int S = memLinks/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
link = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
if (link < memLinks) {
// inside //outside
aq = MassFractionIn; ap = MassFractionOut;
iq = membrane[2*link]; ip = membrane[2*link+1];
nq = iq%Np; np = ip%Np;
nqm = Map[nq]; npm = Map[np]; // strided layout
/* membrane potential for this link */
membranePotential = Psi[nqm] - Psi[npm];
if (membranePotential > Threshold){
aq = ThresholdMassFractionIn; ap = ThresholdMassFractionOut;
}
/* Save the mass transfer coefficients */
coef[2*link] = aq; coef[2*link+1] = ap;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(
const int Cqx, const int Cqy, int const Cqz,
int *Map, double *Distance, double *Psi, double Threshold,
double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count,
const int N, const int Nx, const int Ny, const int Nz) {
//....................................................................................
// Unack distribution from the recv buffer
// Distribution q matche Cqx, Cqy, Cqz
// swap rule means that the distributions in recvbuf are OPPOSITE of q
// dist may be even or odd distributions stored by stream layout
//....................................................................................
int n, idx, link, nqm, npm, i, j, k;
double distanceLocal, distanceNonlocal;
double psiLocal, psiNonlocal, membranePotential;
double ap,aq; // coefficient
/* second enforce custom rule for membrane links */
int S = (count-nlinks)/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
link = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + nlinks;
if (link < count) {
// get the index for the recv list (deal with reordering of links)
idx = d3q7_linkList[link]; // THINK start NEEDS TO BE HERE
// get the distribution index
n = d3q7_recvlist[start+idx];
// get the index in strided layout
nqm = Map[n];
distanceLocal = Distance[nqm];
psiLocal = Psi[nqm];
// Get the 3-D indices from the send process
k = nqm/(Nx*Ny); j = (nqm-Nx*Ny*k)/Nx; i = nqm-Nx*Ny*k-Nx*j;
// Streaming link the non-local distribution
i -= Cqx; j -= Cqy; k -= Cqz;
npm = k*Nx*Ny + j*Nx + i;
distanceNonlocal = Distance[npm];
psiNonlocal = Psi[npm];
membranePotential = psiLocal - psiNonlocal;
aq = MassFractionIn;
ap = MassFractionOut;
/* link is inside membrane */
if (distanceLocal > 0.0){
if (membranePotential < Threshold*(-1.0)){
ap = MassFractionIn;
aq = MassFractionOut;
}
else {
ap = ThresholdMassFractionIn;
aq = ThresholdMassFractionOut;
}
}
else if (membranePotential > Threshold){
aq = ThresholdMassFractionIn;
ap = ThresholdMassFractionOut;
}
// update link based on mass transfer coefficients
coef[2*(link-nlinks)] = aq;
coef[2*(link-nlinks)+1] = ap;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_Membrane_Unpack(int q,
int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count,
double *recvbuf, double *dist, int N, double *coef) {
//....................................................................................
// Unack distribution from the recv buffer
// Distribution q matche Cqx, Cqy, Cqz
// swap rule means that the distributions in recvbuf are OPPOSITE of q
// dist may be even or odd distributions stored by stream layout
//....................................................................................
int n, idx, link;
double fq,fp,fqq,ap,aq; // coefficient
/* second enforce custom rule for membrane links */
int S = count/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
link = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
/* First unpack the regular links */
if (link < nlinks) {
// get the index for the recv list (deal with reordering of links)
idx = d3q7_linkList[link];
// get the distribution index
n = d3q7_recvlist[start+idx];
if (!(n < 0)){
fp = recvbuf[start + idx];
dist[q * N + n] = fp;
}
}
else if (link < count){
/* second enforce custom rule for membrane links */
// get the index for the recv list (deal with reordering of links)
idx = d3q7_linkList[link];
// get the distribution index
n = d3q7_recvlist[start+idx];
// update link based on mass transfer coefficients
if (!(n < 0)){
aq = coef[2*(link-nlinks)];
ap = coef[2*(link-nlinks)+1];
fq = dist[q * N + n];
fp = recvbuf[start + idx];
fqq = (1-aq)*fq+ap*fp;
dist[q * N + n] = fqq;
}
}
}
}
__global__ void dvc_ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef,
double *dist, double *Den, int memLinks, int Np){
int link,iq,ip,nq,np;
double aq, ap, fq, fp, fqq, fpp, Cq, Cp;
int S = memLinks/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
link = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
if (link < memLinks){
// inside //outside
aq = coef[2*link]; ap = coef[2*link+1];
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;
Cq = Den[nq]; Cp = Den[np];
Cq += fqq - fq; Cp += fpp - fp;
Den[nq] = Cq; Den[np] = Cp;
dist[iq] = fqq; dist[ip] = fpp;
}
}
}
__global__ void dvc_ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList, double *dist, double *Den, int start, int finish, int Np){
int n,nread;
@ -107,6 +280,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
double X,Y,Z,factor_x,factor_y,factor_z;
int nr1,nr2,nr3,nr4,nr5,nr6;
int S = Np/NBLOCKS/NTHREADS + 1;
@ -161,34 +335,47 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, doub
FluxElectrical[n+0*Np] = uEPx*Ci;
FluxElectrical[n+1*Np] = uEPy*Ci;
FluxElectrical[n+2*Np] = uEPz*Ci;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);
Z = 4.0 * (uz + uEPz);
factor_x = X / sqrt(1 + X*X);
factor_y = Y / sqrt(1 + Y*Y);
factor_z = Z / sqrt(1 + Z*Z);
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
//dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1
dist[nr2] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx));
//dist[nr2] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr2] =
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
//f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
// q=2
dist[nr1] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx));
//dist[nr1] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr1] =
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
//f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
// q = 3
dist[nr4] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy));
//dist[nr4] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr4] =
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y );
//f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
// q = 4
dist[nr3] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy));
//dist[nr3] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr3] =
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
//f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
// q = 5
dist[nr6] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz));
//dist[nr6] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr6] =
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
//f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
// q = 6
dist[nr5] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz));
//dist[nr5] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[nr5] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
}
}
}
@ -202,6 +389,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *F
double Ex,Ey,Ez;//electrical field
double flux_diffusive_x,flux_diffusive_y,flux_diffusive_z;
double f0,f1,f2,f3,f4,f5,f6;
double X,Y,Z,factor_x,factor_y,factor_z;
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
@ -243,33 +431,46 @@ __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;
/* use logistic function to prevent negative distributions*/
X = 4.0 * (ux + uEPx);
Y = 4.0 * (uy + uEPy);
Z = 4.0 * (uz + uEPz);
factor_x = X / sqrt(1 + X*X);
factor_y = Y / sqrt(1 + Y*Y);
factor_z = Z / sqrt(1 + Z*Z);
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
//dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci*(1.0 - 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[n] = f0 * (1.0 - rlx) + rlx * 0.25 * Ci;
// q = 1
dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx));
//dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[1 * Np + n] =
f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_x);
//f1 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (ux + uEPx));
// q=2
dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx));
//dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx)+8.0*(ux+uEPx)*(ux+uEPx)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[2 * Np + n] =
f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_x);
//f2 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (ux + uEPx));
// q = 3
dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy));
//dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[3 * Np + n] =
f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_y);
//f3 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uy + uEPy));
// q = 4
dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy));
//dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy)+8.0*(uy+uEPy)*(uy+uEPy)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[4 * Np + n] =
f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_y);
//f4 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uy + uEPy));
// q = 5
dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz));
//dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[5 * Np + n] =
f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + factor_z);
//f5 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 + 4.0 * (uz + uEPz));
// q = 6
dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz));
//dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz)+8.0*(uz+uEPz)*(uz+uEPz)- 2.0*((ux+uEPx)*(ux+uEPx) + (uy+uEPy)*(uy+uEPy) + (uz+uEPz)*(uz+uEPz)));
dist[6 * Np + n] =
f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - factor_z);
//f6 * (1.0 - rlx) + rlx * 0.125 * Ci * (1.0 - 4.0 * (uz + uEPz));
}
}
}
@ -420,3 +621,60 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity
}
//cudaProfilerStop();
}
extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef(int *membrane, int *Map, double *Distance, double *Psi, double *coef,
double Threshold, double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int memLinks, int Nx, int Ny, int Nz, int Np){
dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef<<<NBLOCKS,NTHREADS >>>(membrane, Map, Distance, Psi, coef,
Threshold, MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut,
memLinks, Nx, Ny, Nz, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef: %s \n",hipGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo(
const int Cqx, const int Cqy, int const Cqz,
int *Map, double *Distance, double *Psi, double Threshold,
double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn, double ThresholdMassFractionOut,
int *d3q7_recvlist, int *d3q7_linkList, double *coef, int start, int nlinks, int count,
const int N, const int Nx, const int Ny, const int Nz) {
dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo<<<NBLOCKS,NTHREADS >>>(
Cqx, Cqy, Cqz, Map, Distance, Psi, Threshold,
MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut,
d3q7_recvlist, d3q7_linkList, coef, start, nlinks, count, N, Nx, Ny, Nz);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_AssignLinkCoef_halo: %s \n",hipGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_Membrane_Unpack(int q,
int *d3q7_recvlist, int *d3q7_linkList, int start, int nlinks, int count,
double *recvbuf, double *dist, int N, double *coef) {
dvc_ScaLBL_D3Q7_Membrane_Unpack<<<NBLOCKS,NTHREADS >>>(q, d3q7_recvlist, d3q7_linkList, start, nlinks, count,
recvbuf, dist, N, coef) ;
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_Unpack: %s \n",hipGetErrorString(err));
}
}
extern "C" void ScaLBL_D3Q7_Membrane_IonTransport(int *membrane, double *coef,
double *dist, double *Den, int memLinks, int Np){
dvc_ScaLBL_D3Q7_Membrane_IonTransport<<<NBLOCKS,NTHREADS >>>(membrane, coef, dist, Den, memLinks, Np);
hipError_t err = hipGetLastError();
if (hipSuccess != err){
printf("CUDA error in dvc_ScaLBL_D3Q7_Membrane_IonTransport: %s \n",hipGetErrorString(err));
}
}

4
models/ColorModel.cpp
View File

@ -387,6 +387,10 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase) {
AFFINITY, volume_fraction);
}
}
// clean up
delete [] label_count;
delete [] label_count_global;
}
void ScaLBL_ColorModel::Create() {

319
models/IonModel.cpp
View File

@ -583,6 +583,72 @@ void ScaLBL_IonModel::SetDomain() {
nprocz = Dm->nprocz();
}
void ScaLBL_IonModel::SetMembrane() {
size_t NLABELS = 0;
membrane_db = db->getDatabase("Membrane");
/* set distance based on labels inside the membrane (all other labels will be outside) */
auto MembraneLabels = membrane_db->getVector<int>("MembraneLabels");
IonMembrane = std::shared_ptr<Membrane>(new Membrane(Dm, NeighborList, Np));
signed char LABEL = 0;
double *label_count;
double *label_count_global;
Array<char> membrane_id(Nx,Ny,Nz);
label_count = new double[NLABELS];
label_count_global = new double[NLABELS];
// Assign the labels
for (size_t idx = 0; idx < NLABELS; idx++)
label_count[idx] = 0;
/* set the distance to the membrane */
MembraneDistance.resize(Nx, Ny, Nz);
MembraneDistance.fill(0);
for (int k = 0; k < Nz; k++) {
for (int j = 0; j < Ny; j++) {
for (int i = 0; i < Nx; i++) {
membrane_id(i,j,k) = 1; // default value
LABEL = Dm->id[k*Nx*Ny + j*Nx + i];
for (size_t m=0; m<MembraneLabels.size(); m++){
if (LABEL == MembraneLabels[m]) {
label_count[m] += 1.0;
membrane_id(i,j,k) = 0; // inside
m = MembraneLabels.size(); //exit loop
}
}
}
}
}
for (size_t m=0; m<MembraneLabels.size(); m++){
label_count_global[m] = Dm->Comm.sumReduce(label_count[m]);
}
if (rank == 0) {
printf(" Membrane labels: %lu \n", MembraneLabels.size());
for (size_t m=0; m<MembraneLabels.size(); m++){
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);
}
}
/* signed distance to the membrane ( - inside / + outside) */
for (int k = 0; k < Nz; k++) {
for (int j = 0; j < Ny; j++) {
for (int i = 0; i < Nx; i++) {
MembraneDistance(i, j, k) = 2.0 * double(membrane_id(i, j, k)) - 1.0;
}
}
}
CalcDist(MembraneDistance, membrane_id, *Dm);
/* create the membrane data structure */
MembraneCount = IonMembrane->Create(Dm, MembraneDistance, Map);
// clean up
delete [] label_count;
delete [] label_count_global;
}
void ScaLBL_IonModel::ReadInput() {
sprintf(LocalRankString, "%05d", Dm->rank());
@ -778,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,
@ -794,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();
@ -825,7 +923,7 @@ void ScaLBL_IonModel::Initialize() {
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) {
if (File_ion.size() == 2*number_ion_species) {
double *Ci_host;
Ci_host = new double[number_ion_species * Np];
for (size_t ic = 0; ic < number_ion_species; ic++) {
@ -1181,6 +1279,225 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField) {
//if (rank==0) printf("********************************************************\n");
}
void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, double *Psi) {
//Input parameter:
//1. Velocity is from StokesModel
//2. ElectricField is from Poisson model
//LB-related parameter
vector<double> rlx;
for (size_t ic = 0; ic < tau.size(); ic++) {
rlx.push_back(1.0 / tau[ic]);
}
//.......create and start timer............
//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]) {
//************************************************************************/
// *************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
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);
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, 0,
ScaLBL_Comm->LastExterior(), Np);
if (BoundaryConditionSolid == 1) {
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic * Np * 7], IonSolid);
}
ScaLBL_Comm->Barrier();
comm.barrier();
// *************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
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);
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);
if (BoundaryConditionSolid == 1) {
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic * Np * 7], IonSolid);
}
ScaLBL_Comm->Barrier();
comm.barrier();
}
}
//Compute charge density for Poisson equation
for (size_t ic = 0; ic < number_ion_species; ic++) {
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic,
ScaLBL_Comm->FirstInterior(),
ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0,
ScaLBL_Comm->LastExterior(), 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");
}
void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration,
const size_t ic) {
//This function wirte out the data in a normal layout (by aggregating all decomposed domains)

11
models/IonModel.h
View File

@ -1,7 +1,6 @@
/*
* Ion transporte LB Model
*/
#ifndef ScaLBL_IonModel_INC
#define ScaLBL_IonModel_INC
@ -16,6 +15,7 @@
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/Membrane.h"
#include "common/MPI.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
@ -30,10 +30,12 @@ public:
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
void SetDomain();
void SetMembrane();
void ReadInput();
void Create();
void Initialize();
void Run(double *Velocity, double *ElectricField);
void RunMembrane(double *Velocity, double *ElectricField, double *Psi);
void getIonConcentration(DoubleArray &IonConcentration, const size_t ic);
void getIonConcentration_debug(int timestep);
void getIonFluxDiffusive(DoubleArray &IonFlux_x, DoubleArray &IonFlux_y,
@ -88,6 +90,7 @@ public:
IntArray Map;
DoubleArray Distance;
int *NeighborList;
int *dvcMap;
double *fq;
double *Ci;
double *ChargeDensity;
@ -97,6 +100,12 @@ public:
double *FluxDiffusive;
double *FluxAdvective;
double *FluxElectrical;
/* these support membrane capabilities */
std::shared_ptr<Database> membrane_db;
std::shared_ptr<Membrane> IonMembrane;
DoubleArray MembraneDistance;
int MembraneCount; // number of links the cross the membrane
private:
Utilities::MPI comm;

2
tests/CMakeLists.txt
View File

@ -6,6 +6,7 @@ ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
ADD_LBPM_EXECUTABLE( lbpm_greyscale_simulator )
ADD_LBPM_EXECUTABLE( lbpm_greyscaleColor_simulator )
ADD_LBPM_EXECUTABLE( lbpm_electrokinetic_SingleFluid_simulator )
ADD_LBPM_EXECUTABLE( lbpm_cell_simulator )
ADD_LBPM_EXECUTABLE( lbpm_freelee_simulator )
ADD_LBPM_EXECUTABLE( lbpm_freelee_SingleFluidBGK_simulator )
ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator )
@ -58,6 +59,7 @@ ADD_LBPM_TEST( TestTopo3D )
ADD_LBPM_TEST( TestFluxBC )
ADD_LBPM_TEST( TestFlowAdaptor )
ADD_LBPM_TEST( TestMap )
ADD_LBPM_TEST( TestMembrane )
#ADD_LBPM_TEST( TestMRT )
#ADD_LBPM_TEST( TestColorGrad )
ADD_LBPM_TEST( TestWideHalo )

256
tests/TestMembrane.cpp Normal file
View File

@ -0,0 +1,256 @@
//*************************************************************************
// Lattice Boltzmann Simulator for Single Phase Flow in Porous Media
// James E. McCLure
//*************************************************************************
#include <stdio.h>
#include <iostream>
#include <fstream>
#include "common/MPI.h"
#include "common/Membrane.h"
#include "common/ScaLBL.h"
using namespace std;
std::shared_ptr<Database> loadInputs( int nprocs )
{
//auto db = std::make_shared<Database>( "Domain.in" );
auto db = std::make_shared<Database>();
db->putScalar<int>( "BC", 0 );
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 32, 32, 32 } );
db->putScalar<int>( "nspheres", 1 );
db->putVector<double>( "L", { 1, 1, 1 } );
return db;
}
//***************************************************************************************
int main(int argc, char **argv)
{
// Initialize MPI
Utilities::startup( argc, argv );
Utilities::MPI comm( MPI_COMM_WORLD );
int check=0;
{
int i,j,k,n;
int rank = comm.getRank();
if (rank == 0){
printf("********************************************************\n");
printf("Running unit test: TestMembrane \n");
printf("********************************************************\n");
}
// Load inputs
auto db = loadInputs( comm.getSize() );
int Nx = db->getVector<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
auto Dm = std::make_shared<Domain>(db,comm);
Nx += 2;
Ny += 2;
Nz += 2;
int N = Nx*Ny*Nz;
//.......................................................................
int Np = 0;
double distance,radius;
DoubleArray Distance(Nx,Ny,Nz);
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Dm->id[n] = 1;
radius = double(Nx)/4;
distance = sqrt(double((i-0.5*Nx)*(i-0.5*Nx)+ (j-0.5*Ny)*(j-0.5*Ny)+ (k-0.5*Nz)*(k-0.5*Nz)))-radius;
if (distance < 0.0 ){
Dm->id[n] = 1;
}
Distance(i,j,k) = distance;
Np++;
}
}
}
Dm->CommInit();
// Create a communicator for the device (will use optimized layout)
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm(new ScaLBL_Communicator(Dm));
//Create a second communicator based on the regular data layout
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular(new ScaLBL_Communicator(Dm));
if (rank==0){
printf("Total domain size = %i \n",N);
printf("Reduced domain size = %i \n",Np);
}
// LBM variables
if (rank==0) printf ("Set up the neighborlist \n");
int Npad=Np+32;
int neighborSize=18*Npad*sizeof(int);
int *neighborList;
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
//......................device distributions.................................
int *NeighborList;
int *dvcMap;
//...........................................................................
ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Npad);
ScaLBL_CopyToDevice(NeighborList, neighborList, 18*Np*sizeof(int));
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id.data(),Np,1);
comm.barrier();
double *dist;
dist = new double [19*Np];
// Check the neighborlist
printf("Check neighborlist: exterior %i, first interior %i last interior %i \n",ScaLBL_Comm->LastExterior(),ScaLBL_Comm->FirstInterior(),ScaLBL_Comm->LastInterior());
for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
for (int q=0; q<18; q++){
int nn = neighborList[q*Np+idx]%Np;
if (nn>Np) printf("neighborlist error (exterior) at q=%i, idx=%i \n",q,idx);
dist[q*Np + idx] = 0.0;
}
}
for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
for (int q=0; q<18; q++){
int nn = neighborList[q*Np+idx]%Np;
if (nn>Np) printf("neighborlist error (exterior) at q=%i, idx=%i \n",q,idx);
dist[q*Np + idx] = 0.0;
}
}
/* create a membrane data structure */
Membrane M(Dm, NeighborList, Np);
int MembraneCount = M.Create(Dm, Distance, Map);
if (rank==0) printf (" Number of membrane links: %i \n", MembraneCount);
/* create a tagged array to show where the mebrane is*/
double *MembraneLinks;
MembraneLinks = new double [Nx*Ny*Nz];
for (int n=0; n<Nx*Ny*Nz; n++) {
MembraneLinks[n] = 0.0;
}
for (int mlink=0; mlink<MembraneCount; mlink++){
int iq = M.membraneLinks[2*mlink];
int jq = M.membraneLinks[2*mlink+1];
dist[iq] = -1.0; // set these distributions to non-zero
dist[jq] = 1.0;
}
for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){
int idx = Map(i,j,k);
double sum = 0.0;
for (int q=0; q<19; q++){
sum += dist[q*Np + idx];
}
int n = k*Nx*Ny + j*Nx + i;
MembraneLinks[n] = sum;
if (sum > 0.f){
Dm->id[n] = 127;
}
if (sum < 0.f){
Dm->id[n] = 64;
}
}
}
}
if (argc > 1)
Dm->AggregateLabels("membrane.raw");
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("Setting up device map and neighbor list \n");
int *TmpMap;
TmpMap=new int[Np*sizeof(int)];
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
if (!(idx < 0))
TmpMap[idx] = k*Nx*Ny+j*Nx+i;
}
}
}
ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
ScaLBL_DeviceBarrier();
// Create a dummy distribution data structure
double *fq_host;
fq_host = new double[19*Np];
if (rank==0) printf ("Setting up Np=%i distributions \n",Np);
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
if (!(idx<0)){
for (int q=0; q<19; q++){
fq_host[q*Np+idx]=(k*Nx*Ny+j*Nx+i)+0.01*q;
}
}
}
}
}
/* Run dummy communications */
double * gq, *fq;
ScaLBL_AllocateDeviceMemory((void **) &gq, sizeof(double)*7*Np);
ScaLBL_AllocateDeviceMemory((void **) &fq, sizeof(double)*7*Np);
/*initialize fq from host data */
ScaLBL_CopyToDevice(fq, fq_host, sizeof(double)*7*Np);
M.SendD3Q7AA(&fq[0]);
M.RecvD3Q7AA(&gq[0]);
// this has only the communicated values
//ScaLBL_CopyToHost(fq_host, gq, sizeof(double)*7*Np);
if (rank==0) printf ("Sum result \n");
double *Ci;
Ci = new double [Np];
ScaLBL_D3Q7_AAeven_IonConcentration(&gq[0 * Np * 7], &Ci[0 * Np],
0, ScaLBL_Comm->LastExterior(),
Np);
DoubleArray Result(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map, Ci, Result);
/* for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
double sum = 0.0;
if (!(idx<0)){
for (int q=1; q<3; q++){
sum += fq_host[q*Np+idx];
}
Result[k*Nx*Ny+j*Nx+i] = sum;
}
}
}
}
*/
FILE *OUTFILE;
OUTFILE = fopen("D3Q7.raw","wb");
fwrite(Result.data(),8,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
FILE *MAPFILE;
MAPFILE = fopen("Map.raw","wb");
fwrite(Map.data(),4,Nx*Ny*Nz,MAPFILE);
fclose(MAPFILE);
delete [] TmpMap;
delete [] fq_host;
}
Utilities::shutdown();
return check;
}

131
tests/lbpm_cell_simulator.cpp Normal file
View File

@ -0,0 +1,131 @@
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include <math.h>
#include "models/IonModel.h"
#include "models/StokesModel.h"
#include "models/PoissonSolver.h"
#include "models/MultiPhysController.h"
#include "common/Utilities.h"
#include "analysis/ElectroChemistry.h"
using namespace std;
//***************************************************************************
// Test lattice-Boltzmann Ion Model coupled with Poisson equation
//***************************************************************************
int main(int argc, char **argv)
{
// Initialize MPI and error handlers
Utilities::startup( argc, argv );
Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank();
int nprocs = comm.getSize();
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
if (rank == 0){
printf("********************************************************\n");
printf("Running LBPM electrokinetic single-fluid solver \n");
printf("********************************************************\n");
}
// Initialize compute device
int device=ScaLBL_SetDevice(rank);
NULL_USE( device );
ScaLBL_DeviceBarrier();
comm.barrier();
PROFILE_ENABLE(1);
//PROFILE_ENABLE_TRACE();
//PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START("Main");
Utilities::setErrorHandlers();
auto filename = argv[1];
ScaLBL_StokesModel StokesModel(rank,nprocs,comm);
ScaLBL_IonModel IonModel(rank,nprocs,comm);
ScaLBL_Poisson PoissonSolver(rank,nprocs,comm);
ScaLBL_Multiphys_Controller Study(rank,nprocs,comm);//multiphysics controller coordinating multi-model coupling
// Load controller information
Study.ReadParams(filename);
// Load user input database files for Navier-Stokes and Ion solvers
StokesModel.ReadParams(filename);
IonModel.ReadParams(filename);
// Setup other model specific structures
StokesModel.SetDomain();
StokesModel.ReadInput();
StokesModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables
IonModel.SetDomain();
IonModel.ReadInput();
IonModel.Create();
IonModel.SetMembrane();
// Create analysis object
ElectroChemistryAnalyzer Analysis(IonModel.Dm);
// Get internal iteration number
StokesModel.timestepMax = Study.getStokesNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
StokesModel.Initialize(); // initializing the model will set initial conditions for variables
IonModel.timestepMax = Study.getIonNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
IonModel.Initialize();
// Get maximal time converting factor based on Sotkes and Ion solvers
Study.getTimeConvMax_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
// Initialize LB-Poisson model
PoissonSolver.ReadParams(filename);
PoissonSolver.SetDomain();
PoissonSolver.ReadInput();
PoissonSolver.Create();
PoissonSolver.Initialize(Study.time_conv_max);
int timestep=0;
while (timestep < Study.timestepMax){
timestep++;
PoissonSolver.Run(IonModel.ChargeDensity,timestep);//solve Poisson equtaion to get steady-state electrical potental
StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
IonModel.RunMembrane(StokesModel.Velocity,PoissonSolver.ElectricField,PoissonSolver.Psi); //solve for ion transport with membrane
timestep++;//AA operations
if (timestep%Study.analysis_interval==0){
Analysis.Basic(IonModel,PoissonSolver,StokesModel,timestep);
}
if (timestep%Study.visualization_interval==0){
Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep);
/* PoissonSolver.getElectricPotential(timestep);
PoissonSolver.getElectricField(timestep);
IonModel.getIonConcentration(timestep);
StokesModel.getVelocity(timestep);
*/
}
}
if (rank==0) printf("Save simulation raw data at maximum timestep\n");
Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep);
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
if (rank==0) printf("*************************************************************\n");
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_electrokinetic_SingleFluid_simulator",1);
// ****************************************************
} // Limit scope so variables that contain communicators will free before MPI_Finialize
Utilities::shutdown();
}