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

try at membrane simulator

Browse Source
This commit is contained in:
James McClure 2022-03-18 18:08:44 -04:00
parent c94d4e5194
commit 702eaae1c1
5 changed files with 270 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
common/Membrane.cpp
View File

@ -3,21 +3,21 @@
#include "common/Membrane.h"
#include "analysis/distance.h"
Membrane::Membrane(std::shared_ptr <Domain> Dm, int *initialNeighborList, int Nsites) {
Membrane::Membrane(std::shared_ptr <Domain> Dm, int *dvcNeighborList, int Nsites) {
Np = Nsites;
int * initialNeighborList = new int[18*Np];
ScaLBL_AllocateDeviceMemory((void **)&NeighborList, 18*Np*sizeof(int));
Lock=false; // unlock the communicator
//......................................................................................
// Create a separate copy of the communicator for the device
MPI_COMM_SCALBL = Dm->Comm.dup();
Np = Nsites;
neighborList = new int[18*Np];
/* Copy neighborList */
for (int idx=0; idx<Np; idx++){
for (int q = 0; q<18; q++){
neighborList[q*Np+idx] = initialNeighborList[q*18+idx];
}
}
ScaLBL_CopyToHost(initialNeighborList, dvcNeighborList, 18*Np*sizeof(int));
Dm->Comm.barrier();
ScaLBL_CopyToDevice(NeighborList, initialNeighborList, 18*Np*sizeof(int));
/* Copy communication lists */
//......................................................................................
//Lock=false; // unlock the communicator
@ -275,6 +275,11 @@ Membrane::Membrane(std::shared_ptr <Domain> Dm, int *initialNeighborList, int Ns
Membrane::~Membrane() {
delete [] initialNeighborList;
delete [] membraneLinks;
delete [] membraneTag;
delete [] membraneDist;
ScaLBL_FreeDeviceMemory( coefficient_x );
ScaLBL_FreeDeviceMemory( coefficient_X );
ScaLBL_FreeDeviceMemory( coefficient_y );
@ -282,9 +287,9 @@ Membrane::~Membrane() {
ScaLBL_FreeDeviceMemory( coefficient_z );
ScaLBL_FreeDeviceMemory( coefficient_Z );
ScaLBL_FreeDeviceMemory( dvcNeighborList );
ScaLBL_FreeDeviceMemory( dvcMembraneLinks );
ScaLBL_FreeDeviceMemory( dvcMembraneCoef );
ScaLBL_FreeDeviceMemory( NeighborList );
ScaLBL_FreeDeviceMemory( MembraneLinks );
ScaLBL_FreeDeviceMemory( MembraneCoef );
ScaLBL_FreeDeviceMemory( sendbuf_x );
ScaLBL_FreeDeviceMemory( sendbuf_X );
@ -401,6 +406,15 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
int i,j,k;
int idx, neighbor;
double dist, locdist;
int * neighborList = new int[18*Np];
/* Copy neighborList */
for (int idx=0; idx<Np; idx++){
for (int q = 0; q<18; q++){
neighborList[q*Np+idx] = initialNeighborList[q*Np+idx];
}
}
/* go through the neighborlist structure */
/* count & cut the links */
for (k=1;k<Nz-1;k++){
@ -531,7 +545,7 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
}
}
}
/* allocate memory */
membraneTag = new int [mlink];
membraneLinks = new int [2*mlink];
@ -719,14 +733,13 @@ int Membrane::Create(std::shared_ptr <Domain> Dm, DoubleArray &Distance, IntArra
}
/* Create device copies of data structures */
ScaLBL_AllocateDeviceMemory((void **)&dvcNeighborList, 18*Np*sizeof(int));
ScaLBL_AllocateDeviceMemory((void **)&dvcMembraneLinks, 2*mlink*sizeof(int));
ScaLBL_AllocateDeviceMemory((void **)&dvcMembraneCoef, 2*mlink*sizeof(double));
ScaLBL_AllocateDeviceMemory((void **)&dvcMembraneDistance, 2*mlink*sizeof(double));
ScaLBL_AllocateDeviceMemory((void **)&MembraneLinks, 2*mlink*sizeof(int));
ScaLBL_AllocateDeviceMemory((void **)&MembraneCoef, 2*mlink*sizeof(double));
ScaLBL_AllocateDeviceMemory((void **)&MembraneDistance, 2*mlink*sizeof(double));
ScaLBL_CopyToDevice(dvcNeighborList, neighborList, 18*Np*sizeof(int));
ScaLBL_CopyToDevice(dvcMembraneLinks, membraneLinks, 2*mlink*sizeof(int));
ScaLBL_CopyToDevice(dvcMembraneDistance, membraneDist, 2*mlink*sizeof(double));
ScaLBL_CopyToDevice(NeighborList, neighborList, 18*Np*sizeof(int));
ScaLBL_CopyToDevice(MembraneLinks, membraneLinks, 2*mlink*sizeof(int));
ScaLBL_CopyToDevice(MembraneDistance, membraneDist, 2*mlink*sizeof(double));
/* Re-organize communication based on membrane structure*/
@ -1126,7 +1139,7 @@ void Membrane::SendD3Q7AA(double *dist){
req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z, recvCount_z,rank_z,recvtag);
}
void Membrane::RecvD37AA(double *dist){
void Membrane::RecvD3Q7AA(double *dist){
//...................................................................................
// Wait for completion of D3Q19 communication

14
common/Membrane.h
View File

@ -54,7 +54,10 @@ public:
int Np;
int Nx,Ny,Nz,N;
int *neighborList; // modified neighborlist
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
@ -62,10 +65,9 @@ public:
/*
* Device data structures
*/
int *dvcNeighborList;
double *dvcMembraneLinks;
double *dvcMembraneCoef; // mass transport coefficient for the membrane
double *dvcMembraneDistance;
double *MembraneLinks;
double *MembraneCoef; // mass transport coefficient for the membrane
double *MembraneDistance;
/**
* \brief Create a flow adaptor to operate on the LB model
@ -91,7 +93,7 @@ public:
void SendD3Q19AA(double *dist);
void RecvD3Q19AA(double *dist);
void SendD3Q7AA(double *dist);
void RecvD37AA(double *dist);
void RecvD3Q7AA(double *dist);
void AssignCoefficients(int *Map, double *Psi, double *Distance, std::string method);
//......................................................................................
// Buffers to store data sent and recieved by this MPI process

217
models/IonModel.cpp
View File

@ -1246,6 +1246,223 @@ 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();
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],
&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)

17
tests/TestMembrane.cpp
View File

@ -95,6 +95,14 @@ int main(int argc, char **argv)
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();
@ -120,7 +128,7 @@ int main(int argc, char **argv)
}
/* create a membrane data structure */
Membrane M(Dm, neighborList, Np);
Membrane M(Dm, NeighborList, Np);
int MembraneCount = M.Create(Dm, Distance, Map);
if (rank==0) printf (" Number of membrane links: %i \n", MembraneCount);
@ -159,13 +167,6 @@ int main(int argc, char **argv)
if (argc > 1)
Dm->AggregateLabels("membrane.raw");
//......................device distributions.................................
int *NeighborList;
int *dvcMap;
//...........................................................................
ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Npad);
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("Setting up device map and neighbor list \n");

3
tests/lbpm_cell_simulator.cpp
View File

@ -69,6 +69,7 @@ int main(int argc, char **argv)
IonModel.SetDomain();
IonModel.ReadInput();
IonModel.Create();
IonModel.SetMembrane();
// Create analysis object
ElectroChemistryAnalyzer Analysis(IonModel.Dm);
@ -95,7 +96,7 @@ int main(int argc, char **argv)
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.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
IonModel.RunMembrane(StokesModel.Velocity,PoissonSolver.ElectricField,PoissonSolver.Psi); //solve for ion transport with membrane
timestep++;//AA operations