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

add membrane transport function for d3q7

Browse Source
This commit is contained in:
James E McClure
2022-03-07 15:36:16 -05:00
parent dd9da41d4b
commit f270604e6b
4 changed files with 106 additions and 77 deletions
Download Patch File Download Diff File Expand all files Collapse all files

158
analysis/Membrane.cpp
View File

@@ -522,13 +522,18 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
} }
/* allocate memory */ /* allocate memory */
membraneTag = new int [mlink];
membraneLinks = new int [2*mlink]; membraneLinks = new int [2*mlink];
membraneDist = new double [2*mlink]; membraneDist = new double [2*mlink];
membraneTag = new int [mlink]; membraneCoef = new double [2*mlink];
/* construct the membrane*/ /* construct the membrane*/
/* *
* Sites inside the membrane (negative distance) -- store at 2*mlink
* Sites outside the membrane (positive distance) -- store at 2*mlink+1
*/
mlink = 0; mlink = 0;
int insideMem = 0; int outsideMem = 0; int localSite = 0; int neighborSite = 0;
for (k=1;k<Nz-1;k++){ for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){ for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){ for (i=1;i<Nx-1;i++){
@@ -541,17 +546,17 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
dist=Distance(i+1,j,k); dist=Distance(i+1,j,k);
if (dist*locdist < 0.0){ if (dist*locdist < 0.0){
if (locdist < 0.0){ if (locdist < 0.0){
insideMem = 2*mlink; localSite = 2*mlink;
outsideMem = 2*mlink+1; neighborSite = 2*mlink+1;
} }
else{ else{
insideMem = 2*mlink+1; localSite = 2*mlink+1;
outsideMem = 2*mlink; neighborSite = 2*mlink;
} }
membraneLinks[insideMem] = idx + 1*Np; membraneLinks[localSite] = idx + 1*Np;
membraneLinks[outsideMem] = neighbor + 2*Np; membraneLinks[neighborSite] = neighbor + 2*Np;
membraneDist[insideMem] = locdist; membraneDist[localSite] = locdist;
membraneDist[outsideMem] = dist; membraneDist[neighborSite] = dist;
mlink++; mlink++;
} }
@@ -559,17 +564,17 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
dist=Distance(i,j+1,k); dist=Distance(i,j+1,k);
if (dist*locdist < 0.0){ if (dist*locdist < 0.0){
if (locdist < 0.0){ if (locdist < 0.0){
insideMem = 2*mlink; localSite = 2*mlink;
outsideMem = 2*mlink+1; neighborSite = 2*mlink+1;
} }
else{ else{
insideMem = 2*mlink+1; localSite = 2*mlink+1;
outsideMem = 2*mlink; neighborSite = 2*mlink;
} }
membraneLinks[insideMem] = idx + 3*Np; membraneLinks[localSite] = idx + 3*Np;
membraneLinks[outsideMem] = neighbor + 4*Np; membraneLinks[neighborSite] = neighbor + 4*Np;
membraneDist[insideMem] = locdist; membraneDist[localSite] = locdist;
membraneDist[outsideMem] = dist; membraneDist[neighborSite] = dist;
mlink++; mlink++;
} }
@@ -577,17 +582,17 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
dist=Distance(i,j,k+1); dist=Distance(i,j,k+1);
if (dist*locdist < 0.0){ if (dist*locdist < 0.0){
if (locdist < 0.0){ if (locdist < 0.0){
insideMem = 2*mlink; localSite = 2*mlink;
outsideMem = 2*mlink+1; neighborSite = 2*mlink+1;
} }
else{ else{
insideMem = 2*mlink+1; localSite = 2*mlink+1;
outsideMem = 2*mlink; neighborSite = 2*mlink;
} }
membraneLinks[insideMem] = idx + 5*Np; membraneLinks[localSite] = idx + 5*Np;
membraneLinks[outsideMem] = neighbor + 6*Np; membraneLinks[neighborSite] = neighbor + 6*Np;
membraneDist[insideMem] = locdist; membraneDist[localSite] = locdist;
membraneDist[outsideMem] = dist; membraneDist[neighborSite] = dist;
mlink++; mlink++;
} }
@@ -595,17 +600,17 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
dist=Distance(i+1,j+1,k); dist=Distance(i+1,j+1,k);
if (dist*locdist < 0.0){ if (dist*locdist < 0.0){
if (locdist < 0.0){ if (locdist < 0.0){
insideMem = 2*mlink; localSite = 2*mlink;
outsideMem = 2*mlink+1; neighborSite = 2*mlink+1;
} }
else{ else{
insideMem = 2*mlink+1; localSite = 2*mlink+1;
outsideMem = 2*mlink; neighborSite = 2*mlink;
} }
membraneLinks[insideMem] = idx + 7*Np; membraneLinks[localSite] = idx + 7*Np;
membraneLinks[outsideMem] = neighbor+8*Np; membraneLinks[neighborSite] = neighbor+8*Np;
membraneDist[insideMem] = locdist; membraneDist[localSite] = locdist;
membraneDist[outsideMem] = dist; membraneDist[neighborSite] = dist;
mlink++; mlink++;
} }
@@ -613,17 +618,17 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
dist=Distance(i+1,j-1,k); dist=Distance(i+1,j-1,k);
if (dist*locdist < 0.0){ if (dist*locdist < 0.0){
if (locdist < 0.0){ if (locdist < 0.0){
insideMem = 2*mlink; localSite = 2*mlink;
outsideMem = 2*mlink+1; neighborSite = 2*mlink+1;
} }
else{ else{
insideMem = 2*mlink+1; localSite = 2*mlink+1;
outsideMem = 2*mlink; neighborSite = 2*mlink;
} }
membraneLinks[insideMem] = idx + 9*Np; membraneLinks[localSite] = idx + 9*Np;
membraneLinks[outsideMem] = neighbor + 10*Np; membraneLinks[neighborSite] = neighbor + 10*Np;
membraneDist[insideMem] = locdist; membraneDist[localSite] = locdist;
membraneDist[outsideMem] = dist; membraneDist[neighborSite] = dist;
mlink++; mlink++;
} }
@@ -631,17 +636,17 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
dist=Distance(i+1,j,k+1); dist=Distance(i+1,j,k+1);
if (dist*locdist < 0.0){ if (dist*locdist < 0.0){
if (locdist < 0.0){ if (locdist < 0.0){
insideMem = 2*mlink; localSite = 2*mlink;
outsideMem = 2*mlink+1; neighborSite = 2*mlink+1;
} }
else{ else{
insideMem = 2*mlink+1; localSite = 2*mlink+1;
outsideMem = 2*mlink; neighborSite = 2*mlink;
} }
membraneLinks[insideMem] = idx + 11*Np; membraneLinks[localSite] = idx + 11*Np;
membraneLinks[outsideMem] = neighbor + 12*Np; membraneLinks[neighborSite] = neighbor + 12*Np;
membraneDist[insideMem] = locdist; membraneDist[localSite] = locdist;
membraneDist[outsideMem] = dist; membraneDist[neighborSite] = dist;
mlink++; mlink++;
} }
@@ -649,17 +654,17 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
dist=Distance(i+1,j,k-1); dist=Distance(i+1,j,k-1);
if (dist*locdist < 0.0){ if (dist*locdist < 0.0){
if (locdist < 0.0){ if (locdist < 0.0){
insideMem = 2*mlink; localSite = 2*mlink;
outsideMem = 2*mlink+1; neighborSite = 2*mlink+1;
} }
else{ else{
insideMem = 2*mlink+1; localSite = 2*mlink+1;
outsideMem = 2*mlink; neighborSite = 2*mlink;
} }
membraneLinks[insideMem] = idx + 13*Np; membraneLinks[localSite] = idx + 13*Np;
membraneLinks[outsideMem] = neighbor + 14*Np; membraneLinks[neighborSite] = neighbor + 14*Np;
membraneDist[insideMem] = locdist; membraneDist[localSite] = locdist;
membraneDist[outsideMem] = dist; membraneDist[neighborSite] = dist;
mlink++; mlink++;
} }
@@ -667,17 +672,17 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
dist=Distance(i,j+1,k+1); dist=Distance(i,j+1,k+1);
if (dist*locdist < 0.0){ if (dist*locdist < 0.0){
if (locdist < 0.0){ if (locdist < 0.0){
insideMem = 2*mlink; localSite = 2*mlink;
outsideMem = 2*mlink+1; neighborSite = 2*mlink+1;
} }
else{ else{
insideMem = 2*mlink+1; localSite = 2*mlink+1;
outsideMem = 2*mlink; neighborSite = 2*mlink;
} }
membraneLinks[insideMem] = idx + 15*Np; membraneLinks[localSite] = idx + 15*Np;
membraneLinks[outsideMem] =neighbor + 16*Np; membraneLinks[neighborSite] = neighbor + 16*Np;
membraneDist[insideMem] = locdist; membraneDist[localSite] = locdist;
membraneDist[outsideMem] = dist; membraneDist[neighborSite] = dist;
mlink++; mlink++;
} }
@@ -685,18 +690,17 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
dist=Distance(i,j+1,k-1); dist=Distance(i,j+1,k-1);
if (dist*locdist < 0.0){ if (dist*locdist < 0.0){
if (locdist < 0.0){ if (locdist < 0.0){
insideMem = 2*mlink; localSite = 2*mlink;
outsideMem = 2*mlink+1; neighborSite = 2*mlink+1;
} }
else{ else{
insideMem = 2*mlink+1; localSite = 2*mlink+1;
neighborSite = 2*mlink;
outsideMem = 2*mlink;
} }
membraneLinks[insideMem] = idx + 17*Np; membraneLinks[localSite] = idx + 17*Np;
membraneLinks[outsideMem] = neighbor + 18*Np; membraneLinks[neighborSite] = neighbor + 18*Np;
membraneDist[insideMem] = locdist; membraneDist[localSite] = locdist;
membraneDist[outsideMem] = dist; membraneDist[neighborSite] = dist;
mlink++; mlink++;
} }
} }
@@ -706,10 +710,10 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
/* Re-organize communication based on membrane structure*/ /* Re-organize communication based on membrane structure*/
//...Map recieve list for the X face: q=2,8,10,12,14 ................................. //...Map recieve list for the X face: q=2,8,10,12,14 .................................
linkCount_X[0]= D3Q19_MapRecv(-1,0,0, Dm->recvList("X"),0,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); linkCount_X[0] = D3Q19_MapRecv(-1,0,0, Dm->recvList("X"),0,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map);
linkCount_X[1] = D3Q19_MapRecv(-1,-1,0,Dm->recvList("X"),recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); linkCount_X[1] = D3Q19_MapRecv(-1,-1,0,Dm->recvList("X"),recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map);
linkCount_X[2] = D3Q19_MapRecv(-1,1,0, Dm->recvList("X"),2*recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); linkCount_X[2] = D3Q19_MapRecv(-1,1,0, Dm->recvList("X"),2*recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map);
linkCount_X[3]= D3Q19_MapRecv(-1,0,-1,Dm->recvList("X"),3*recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); linkCount_X[3] = D3Q19_MapRecv(-1,0,-1,Dm->recvList("X"),3*recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map);
linkCount_X[4] = D3Q19_MapRecv(-1,0,1, Dm->recvList("X"),4*recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map); linkCount_X[4] = D3Q19_MapRecv(-1,0,1, Dm->recvList("X"),4*recvCount_X,recvCount_X,dvcRecvDist_X,dvcRecvLinks_X,Distance,Map);
//................................................................................... //...................................................................................
//...Map recieve list for the x face: q=1,7,9,11,13.................................. //...Map recieve list for the x face: q=1,7,9,11,13..................................

1
analysis/Membrane.h
View File

@@ -58,6 +58,7 @@ public:
int *membraneLinks; // D3Q19 links that cross membrane int *membraneLinks; // D3Q19 links that cross membrane
int *membraneTag; // label each link in the membrane int *membraneTag; // label each link in the membrane
double *membraneDist; // distance to membrane for each linked site double *membraneDist; // distance to membrane for each linked site
double *membraneCoef; // mass transport coefficient for the membrane
/** /**
* \brief Create a flow adaptor to operate on the LB model * \brief Create a flow adaptor to operate on the LB model

5
common/ScaLBL.h
View File

@@ -217,6 +217,11 @@ 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); 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);
// GREYSCALE MODEL (Single-component) // GREYSCALE MODEL (Single-component)
extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *Dist, int Np, double Den); extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *Dist, int Np, double Den);

19
cpu/Ion.cpp
View File

@@ -1,5 +1,24 @@
#include <stdio.h> #include <stdio.h>
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, extern "C" void ScaLBL_D3Q7_AAodd_IonConcentration(int *neighborList,
double *dist, double *Den, double *dist, double *Den,
int start, int finish, int start, int finish,