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

it all builds again

Browse Source
This commit is contained in:
James E McClure 2018-05-18 16:07:26 -04:00
parent 64a670b4e1
commit 2e38569a83
7 changed files with 132 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -158,7 +158,7 @@ IF ( NOT ONLY_BUILD_DOCS )
ADD_PACKAGE_SUBDIRECTORY( analysis )
ADD_PACKAGE_SUBDIRECTORY( IO )
ADD_PACKAGE_SUBDIRECTORY( threadpool )
#ADD_PACKAGE_SUBDIRECTORY( models )
ADD_PACKAGE_SUBDIRECTORY( models )
IF ( USE_CUDA )
ADD_PACKAGE_SUBDIRECTORY( gpu )
ELSE()

2
tests/CMakeLists.txt
View File

@ -24,7 +24,7 @@ ADD_LBPM_EXECUTABLE( lbpm_porenetwork_pp )
ADD_LBPM_EXECUTABLE( lbpm_plates_pp )
ADD_LBPM_EXECUTABLE( lbpm_squaretube_pp )
ADD_LBPM_EXECUTABLE( GenerateSphereTest )
ADD_LBPM_EXECUTABLE( ComponentLabel )
#ADD_LBPM_EXECUTABLE( ComponentLabel )
ADD_LBPM_EXECUTABLE( ColorToBinary )
ADD_LBPM_EXECUTABLE( DataAggregator )
ADD_LBPM_EXECUTABLE( BlobAnalysis )

18
tests/TestBlobAnalyze.cpp
View File

@ -56,7 +56,7 @@ inline void WriteBlobs(TwoPhase Averages){
inline void WriteBlobStates(TwoPhase TCAT, double D, double porosity){
int a;
double iVol=1.0/TCAT.Dm.Volume;
double iVol=1.0/TCAT.Dm->Volume;
double PoreVolume;
double nwp_volume,vol_n,pan,pn,pw,pawn,pwn,awn,ans,aws,Jwn,Kwn,lwns,cwns,clwns;
double sw,awnD,awsD,ansD,lwnsDD,JwnD,pc;
@ -216,7 +216,7 @@ int main(int argc, char **argv)
// Set up the domain
int BC=0;
// Get the rank info
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
std::shared_ptr<Domain> Dm(new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
// const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
TwoPhase Averages(Dm);
@ -229,12 +229,12 @@ int main(int argc, char **argv)
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
Dm.id[n] = 1;
Dm->id[n] = 1;
}
}
}
//.......................................................................
Dm.CommInit(comm); // Initialize communications for domains
Dm->CommInit(comm); // Initialize communications for domains
//.......................................................................
// Read in sphere pack (initialize the non-wetting phase as inside of spheres)
//
@ -264,7 +264,7 @@ int main(int argc, char **argv)
MPI_Barrier(comm);
//.......................................................................
SignedDistance(Averages.Phase.data(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
Dm.iproc(),Dm.jproc(),Dm.kproc(),Dm.nprocx(),Dm.nprocy(),Dm.nprocz());
Dm->iproc(),Dm->jproc(),Dm->kproc(),Dm->nprocx(),Dm->nprocy(),Dm->nprocz());
//.......................................................................
// Assign the phase ID field based on the signed distance
//.......................................................................
@ -277,10 +277,10 @@ int main(int argc, char **argv)
Averages.SDs(i,j,k) = 100.0;
Averages.Phase(i,j,k) += 2.0;
if (Averages.Phase(i,j,k) > 0.0){
Dm.id[n] = 2;
Dm->id[n] = 2;
}
else{
Dm.id[n] = 1;
Dm->id[n] = 1;
}
Averages.SDn(i,j,k) = -Averages.Phase(i,j,k);
Averages.Phase(i,j,k) = Averages.SDn(i,j,k);
@ -297,8 +297,8 @@ int main(int argc, char **argv)
if (rank==0) printf("initializing the system \n");
Averages.UpdateSolid();
Dm.CommunicateMeshHalo(Averages.Phase);
Dm.CommunicateMeshHalo(Averages.SDn);
Dm->CommunicateMeshHalo(Averages.Phase);
Dm->CommunicateMeshHalo(Averages.SDn);
Averages.Initialize();
Averages.UpdateMeshValues();

122
tests/TestBubbleDFH.cpp
View File

@ -167,14 +167,14 @@ int main(int argc, char **argv)
pBC=false;
// Full domain used for averaging (do not use mask for analysis)
Domain Dm(domain_db);
for (int i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
std::shared_ptr<Domain> Dm(new Domain(domain_db));
for (int i=0; i<Dm->Nx*Dm->Ny*Dm->Nz; i++) Dm->id[i] = 1;
std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
// TwoPhase Averages(Dm);
Dm.CommInit(comm);
Dm->CommInit(comm);
// Mask that excludes the solid phase
Domain Mask(domain_db);
std::shared_ptr<Domain> Mask(new Domain(domain_db));
MPI_Barrier(comm);
Nx+=2; Ny+=2; Nz += 2;
@ -226,9 +226,9 @@ int main(int argc, char **argv)
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
int n = k*Nx*Ny + j*Nz + i;
int iglobal= i+(Nx-2)*Dm.iproc();
int jglobal= j+(Ny-2)*Dm.jproc();
int kglobal= k+(Nz-2)*Dm.kproc();
int iglobal= i+(Nx-2)*Dm->iproc();
int jglobal= j+(Ny-2)*Dm->jproc();
int kglobal= k+(Nz-2)*Dm->kproc();
// Initialize phase position field for parallel bubble test
if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
+(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)
@ -249,22 +249,22 @@ int main(int argc, char **argv)
//.........................................................
// Initialize communication structures in averaging domain
for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
Mask.CommInit(comm);
for (i=0; i<Mask->Nx*Mask->Ny*Mask->Nz; i++) Mask->id[i] = id[i];
Mask->CommInit(comm);
double *PhaseLabel;
PhaseLabel = new double[N];
//...........................................................................
if (rank==0) printf ("Create ScaLBL_Communicator \n");
// Create a communicator for the device (will use optimized layout)
ScaLBL_Communicator ScaLBL_Comm(Mask);
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm(new ScaLBL_Communicator(Mask));
int Npad=(Np/16 + 2)*16;
if (rank==0) printf ("Set up memory efficient layout Npad=%i \n",Npad);
int *neighborList;
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Mask.id,Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
MPI_Barrier(comm);
//...........................................................................
@ -339,8 +339,8 @@ int main(int argc, char **argv)
Tmp[idx+Np] = value*dy;
Tmp[idx+2*Np] = value*dz;
// initialize fluid phases
if (Mask.id[n] == 1) PhaseLabel[idx] = 1.0;
else if (Mask.id[n] == 2){
if (Mask->id[n] == 1) PhaseLabel[idx] = 1.0;
else if (Mask->id[n] == 2){
PhaseLabel[idx] = -1.0;
count_wet +=1.0;
}
@ -365,13 +365,13 @@ int main(int argc, char **argv)
if (rank==0) printf ("Initializing distributions \n");
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf ("Initializing phase field \n");
ScaLBL_DFH_Init(Phi, Den, Aq, Bq, 0, ScaLBL_Comm.last_interior, Np);
ScaLBL_DFH_Init(Phi, Den, Aq, Bq, 0, ScaLBL_Comm->last_interior, Np);
//.......................................................................
// Once phase has been initialized, map solid to account for 'smeared' interface
//for (i=0; i<N; i++) Averages.SDs(i) -= (1.0);
// Make sure the id match for the two domains
for (i=0; i<N; i++) Dm.id[i] = Mask.id[i];
for (i=0; i<N; i++) Dm->id[i] = Mask->id[i];
//.......................................................................
// Finalize setup for averaging domain
Averages->UpdateSolid();
@ -389,10 +389,10 @@ int main(int argc, char **argv)
//...........................................................................
ScaLBL_DeviceBarrier();
ScaLBL_CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
ScaLBL_Comm.RegularLayout(Map,Pressure,Averages->Press);
ScaLBL_Comm.RegularLayout(Map,&Velocity[0],Averages->Vel_x);
ScaLBL_Comm.RegularLayout(Map,&Velocity[Np],Averages->Vel_y);
ScaLBL_Comm.RegularLayout(Map,&Velocity[2*Np],Averages->Vel_z);
ScaLBL_Comm->RegularLayout(Map,Pressure,Averages->Press);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Averages->Vel_x);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Averages->Vel_y);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Averages->Vel_z);
//...........................................................................
if (rank==0) printf("********************************************************\n");
@ -422,73 +422,73 @@ int main(int argc, char **argv)
timestep++;
// Compute the Phase indicator field
// Read for Aq, Bq happens in this routine (requires communication)
ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, 0, ScaLBL_Comm.next, Np);
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->next, Np);
// compute the gradient
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
ScaLBL_Comm.SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm.next, Np);
ScaLBL_Comm.RecvGrad(Phi,Gradient);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_Comm->SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm->next, Np);
ScaLBL_Comm->RecvGrad(Phi,Gradient);
// Perform the collision operation
ScaLBL_Comm.SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set BCs
if (BoundaryCondition > 0){
ScaLBL_Comm.Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm.Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
if (BoundaryCondition == 3){
ScaLBL_Comm.D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
if (BoundaryCondition == 4){
din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm.next, Np);
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm->next, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// *************EVEN TIMESTEP*************
timestep++;
// Compute the Phase indicator field
ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, 0, ScaLBL_Comm.next, Np);
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, 0, ScaLBL_Comm->next, Np);
// compute the gradient
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
ScaLBL_Comm.SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm.next, Np);
ScaLBL_Comm.RecvGrad(Phi,Gradient);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_Comm->SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm->next, Np);
ScaLBL_Comm->RecvGrad(Phi,Gradient);
// Perform the collision operation
ScaLBL_Comm.SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np);
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
if (BoundaryCondition > 0){
ScaLBL_Comm.Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm.Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
if (BoundaryCondition == 3){
ScaLBL_Comm.D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
else if (BoundaryCondition == 4){
din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm.next, Np);
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm->next, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************
MPI_Barrier(comm);
@ -522,7 +522,7 @@ int main(int argc, char **argv)
// Copy back final phase indicator field and convert to regular layout
DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm.RegularLayout(Map,Phi,PhaseField);
ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField);
FILE *OUTFILE;
sprintf(LocalRankFilename,"Phase.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
@ -533,9 +533,9 @@ int main(int argc, char **argv)
DoubleArray Cy(Nx,Ny,Nz);
DoubleArray Cz(Nx,Ny,Nz);
DoubleArray GradNorm(Nx,Ny,Nz);
ScaLBL_Comm.RegularLayout(Map,&Gradient[0],Cx);
ScaLBL_Comm.RegularLayout(Map,&Gradient[Np],Cy);
ScaLBL_Comm.RegularLayout(Map,&Gradient[2*Np],Cz);
ScaLBL_Comm->RegularLayout(Map,&Gradient[0],Cx);
ScaLBL_Comm->RegularLayout(Map,&Gradient[Np],Cy);
ScaLBL_Comm->RegularLayout(Map,&Gradient[2*Np],Cz);
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
@ -551,8 +551,8 @@ int main(int argc, char **argv)
DoubleArray Rho1(Nx,Ny,Nz);
DoubleArray Rho2(Nx,Ny,Nz);
ScaLBL_Comm.RegularLayout(Map,&Den[0],Rho1);
ScaLBL_Comm.RegularLayout(Map,&Den[Np],Rho2);
ScaLBL_Comm->RegularLayout(Map,&Den[0],Rho1);
ScaLBL_Comm->RegularLayout(Map,&Den[Np],Rho2);
FILE *RFILE1;
sprintf(LocalRankFilename,"Rho1.raw",rank);
RFILE1 = fopen(LocalRankFilename,"wb");

123
tests/TestColorMassBounceback.cpp
View File

@ -27,7 +27,6 @@ int main(int argc, char **argv)
int check;
{
// parallel domain size (# of sub-domains)
int iproc,jproc,kproc;
int i,j,k,n,Npad;
auto filename = argv[1];
auto db = std::make_shared<Database>( filename );
@ -50,7 +49,6 @@ int main(int argc, char **argv)
string FILENAME;
// Color Model parameters
int timestepMax = color_db->getScalar<int>( "timestepMax" );
double tauA = color_db->getScalar<double>( "tauA" );
double tauB = color_db->getScalar<double>( "tauB" );
double rhoA = color_db->getScalar<double>( "rhoA" );
@ -62,32 +60,20 @@ int main(int argc, char **argv)
double beta = color_db->getScalar<double>( "beta" );
bool Restart = color_db->getScalar<bool>( "Restart" );
double din = color_db->getScalar<double>( "din" );
double dout = color_db->getScalar<double>( "dout" );;
double inletA=1.f;
double inletB=0.f;
double outletA=0.f;
double outletB=1.f;
double flux = 10.f;
// Read domain values
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
int BoundaryCondition = domain_db->getScalar<int>( "BC" );
int Nx = size[0];
int Ny = size[1];
int Nz = size[2];
double Lx = L[0];
double Ly = L[1];
double Lz = L[2];
int nprocx = nproc[0];
int nprocy = nproc[1];
int nprocz = nproc[2];
int timestep = 6;
if (BoundaryCondition==4) flux = din*rhoA; // mass flux must adjust for density (see formulation for details
// Get the rank info
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
@ -132,12 +118,9 @@ int main(int argc, char **argv)
pBC=false;
// Full domain used for averaging (do not use mask for analysis)
Domain Dm(domain_db);
for (int i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
Dm.CommInit(comm);
// Mask that excludes the solid phase
Domain Mask(domain_db);
std::shared_ptr<Domain> Dm(new Domain(domain_db));
for (int i=0; i<Dm->Nx*Dm->Ny*Dm->Nz; i++) Dm->id[i] = 1;
Dm->CommInit(comm);
MPI_Barrier(comm);
Nx+=2; Ny+=2; Nz += 2;
@ -265,30 +248,30 @@ int main(int argc, char **argv)
if (rank==0) printf ("Initializing distributions \n");
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf ("Initializing phase field \n");
ScaLBL_DFH_Init(Phi, Den, Aq, Bq, 0, ScaLBL_Comm.LastInterior(), Np);
ScaLBL_DFH_Init(Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastInterior(), Np);
// *************ODD TIMESTEP*************
// Compute the Phase indicator field
// Read for Aq, Bq happens in this routine (requires communication)
// Read for Aq, Bq happens in this routine (requires communication)
ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, 0, ScaLBL_Comm.LastExterior(), Np);
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// compute the gradient
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm.LastExterior(), Np);
ScaLBL_Comm.RecvGrad(Phi,Gradient);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->RecvGrad(Phi,Gradient);
// Perform the collision operation
ScaLBL_Comm.SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm.LastExterior(), Np);
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
@ -331,24 +314,24 @@ int main(int argc, char **argv)
}
// EVEN TIMESTEP
// Compute the Phase indicator field
ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, 0, ScaLBL_Comm.LastExterior(), Np);
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// compute the gradient
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm.LastExterior(), Np);
ScaLBL_Comm.RecvGrad(Phi,Gradient);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->RecvGrad(Phi,Gradient);
// Perform the collision operation
ScaLBL_Comm.SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm.LastExterior(), Np);
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
//************************************************************************
@ -407,31 +390,31 @@ int main(int argc, char **argv)
if (rank==0) printf ("Initializing distributions \n");
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf ("Initializing phase field \n");
ScaLBL_DFH_Init(Phi, Den, Aq, Bq, 0, ScaLBL_Comm.LastInterior(), Np);
ScaLBL_DFH_Init(Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastInterior(), Np);
// *************ODD TIMESTEP*************
// Compute the Phase indicator field
// Read for Aq, Bq happens in this routine (requires communication)
// Read for Aq, Bq happens in this routine (requires communication)
ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, 0, ScaLBL_Comm.LastExterior(), Np);
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAodd_DFH(NeighborList, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// compute the gradient
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm.LastExterior(), Np);
ScaLBL_Comm.RecvGrad(Phi,Gradient);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->RecvGrad(Phi,Gradient);
// Perform the collision operation
ScaLBL_Comm.SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAodd_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm.LastExterior(), Np);
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
@ -475,24 +458,24 @@ int main(int argc, char **argv)
// EVEN TIMESTEP
// Compute the Phase indicator field
ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, 0, ScaLBL_Comm.LastExterior(), Np);
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAeven_DFH(Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// compute the gradient
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm.LastExterior(), Np);
ScaLBL_Comm.RecvGrad(Phi,Gradient);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->SendHalo(Phi);
ScaLBL_D3Q19_Gradient_DFH(NeighborList, Phi, Gradient, SolidPotential, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->RecvGrad(Phi,Gradient);
// Perform the collision operation
ScaLBL_Comm.SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm.FirstInterior(), ScaLBL_Comm.LastInterior(), Np);
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
alpha, beta, Fx, Fy, Fz, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAeven_DFH(NeighborList, fq, Aq, Bq, Den, Phi, Gradient, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm.LastExterior(), Np);
alpha, beta, Fx, Fy, Fz, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
//************************************************************************

20
tests/TestMap.cpp
View File

@ -23,24 +23,18 @@ int main(int argc, char **argv)
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
int check;
int check=0;
{
// parallel domain size (# of sub-domains)
int nprocx,nprocy,nprocz;
int iproc,jproc,kproc;
if (rank == 0){
printf("********************************************************\n");
printf("Running Color Model: TestColor \n");
printf("Running Color Model: TestMap \n");
printf("********************************************************\n");
}
// BGK Model parameters
string FILENAME;
unsigned int nBlocks, nthreads;
int timestepMax, interval;
double Fx,Fy,Fz,tol;
// Domain variables
double Lx,Ly,Lz;
int nspheres;
@ -176,15 +170,15 @@ int main(int argc, char **argv)
MPI_Barrier(comm);
// 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++){
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);
}
}
for (int idx=ScaLBL_Comm.FirstInterior(); idx<ScaLBL_Comm.LastInterior(); idx++){
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);
@ -193,12 +187,8 @@ int main(int argc, char **argv)
}
//......................device distributions.................................
int dist_mem_size = Np*sizeof(double);
if (rank==0) printf ("Allocating distributions \n");
int *NeighborList;
int *dvcMap;
//...........................................................................
ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Npad);

7
tests/lbpm_segmented_pp.cpp
View File

@ -44,7 +44,6 @@ double ReadFromBlock( char *ID, int iproc, int jproc, int kproc, int Nx, int Ny,
char sx[2];
char sy[2];
char sz[2];
char tmpstr[10];
// array to store ids read from block
char *id;
@ -62,9 +61,8 @@ double ReadFromBlock( char *ID, int iproc, int jproc, int kproc, int Nx, int Ny,
//fflush(stdout);
// Read the file
size_t readID;
FILE *IDFILE = fopen(LocalRankFilename,"rb");
readID=fread(id,1,1024*1024*1024,IDFILE);
fread(id,1,1024*1024*1024,IDFILE);
fclose(IDFILE);
//printf("Loading data ... \n");
@ -128,7 +126,6 @@ int main(int argc, char **argv)
//.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int Nx,Ny,Nz;
int i,j,k,n;
int BC=0;
// char fluidValue,solidValue;
@ -254,7 +251,7 @@ int main(int argc, char **argv)
MeanFilter(Averages->SDs);
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
CalcDist(Averages.SDs,id,Dm);
CalcDist(Averages->SDs,id,*Dm);
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");