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

down the rabbit hole

Browse Source
This commit is contained in:
James E McClure
2018-05-16 20:54:37 -04:00
parent 0b8c5f1053
commit 5060f59323
8 changed files with 359 additions and 346 deletions
Download Patch File Download Diff File Expand all files Collapse all files

242
analysis/TwoPhase.cpp
View File

@@ -54,7 +54,7 @@
#define PI 3.14159265359
// Constructor
TwoPhase::TwoPhase(Domain &dm):
TwoPhase::TwoPhase(Domain *dm):
n_nw_pts(0), n_ns_pts(0), n_ws_pts(0), n_nws_pts(0), n_local_sol_pts(0), n_local_nws_pts(0),
n_nw_tris(0), n_ns_tris(0), n_ws_tris(0), n_nws_seg(0), n_local_sol_tris(0),
nc(0), kstart(0), kfinish(0), fluid_isovalue(0), solid_isovalue(0), Volume(0),
@@ -69,8 +69,8 @@ TwoPhase::TwoPhase(Domain &dm):
trRwn(0), trRwn_global(0), nwp_volume_global(0), wp_volume_global(0),
As_global(0), wwndnw_global(0), wwnsdnwn_global(0), Jwnwwndnw_global(0), dEs(0), dAwn(0), dAns(0)
{
Nx=dm.Nx; Ny=dm.Ny; Nz=dm.Nz;
Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm.nprocx()*Dm.nprocy()*Dm.nprocz()*1.0;
Nx=dm->Nx; Ny=dm->Ny; Nz=dm->Nz;
Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm->nprocx()*Dm->nprocy()*Dm->nprocz()*1.0;
TempID = new char[Nx*Ny*Nz];
@@ -135,7 +135,7 @@ TwoPhase::TwoPhase(Domain &dm):
Gns_global.resize(6);
Gws_global.resize(6);
//.........................................
if (Dm.rank()==0){
if (Dm->rank()==0){
TIMELOG = fopen("timelog.tcat","a+");
if (fseek(TIMELOG,0,SEEK_SET) == fseek(TIMELOG,0,SEEK_CUR))
{
@@ -165,7 +165,7 @@ TwoPhase::TwoPhase(Domain &dm):
}
else{
char LocalRankString[8];
sprintf(LocalRankString,"%05d",Dm.rank());
sprintf(LocalRankString,"%05d",Dm->rank());
char LocalRankFilename[40];
sprintf(LocalRankFilename,"%s%s","timelog.tcat.",LocalRankString);
TIMELOG = fopen(LocalRankFilename,"a+");
@@ -259,7 +259,7 @@ void TwoPhase::ComputeDelPhi()
int i,j,k;
double fx,fy,fz;
Dm.CommunicateMeshHalo(Phase);
Dm->CommunicateMeshHalo(Phase);
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
@@ -307,8 +307,8 @@ void TwoPhase::SetupCubes(Domain &Dm)
int i,j,k;
kstart = 1;
kfinish = Nz-1;
if (Dm.BoundaryCondition !=0 && Dm.kproc==0) kstart = 4;
if (Dm.BoundaryCondition !=0 && Dm.kproc==Dm.nprocz-1) kfinish = Nz-4;
if (Dm->BoundaryCondition !=0 && Dm->kproc==0) kstart = 4;
if (Dm->BoundaryCondition !=0 && Dm->kproc==Dm->nprocz-1) kfinish = Nz-4;
nc=0;
for (k=kstart; k<kfinish; k++){
for (j=1; j<Ny-1; j++){
@@ -327,17 +327,17 @@ void TwoPhase::SetupCubes(Domain &Dm)
void TwoPhase::UpdateSolid()
{
Dm.CommunicateMeshHalo(SDs);
Dm->CommunicateMeshHalo(SDs);
//...........................................................................
// Gradient of the Signed Distance function
//...........................................................................
pmmc_MeshGradient(SDs,SDs_x,SDs_y,SDs_z,Nx,Ny,Nz);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_x);
Dm->CommunicateMeshHalo(SDs_x);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_y);
Dm->CommunicateMeshHalo(SDs_y);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_z);
Dm->CommunicateMeshHalo(SDs_z);
//...........................................................................
}
@@ -346,27 +346,27 @@ void TwoPhase::UpdateMeshValues()
{
int i,j,k,n;
//...........................................................................
Dm.CommunicateMeshHalo(SDn);
Dm->CommunicateMeshHalo(SDn);
//...........................................................................
// Compute the gradients of the phase indicator and signed distance fields
pmmc_MeshGradient(SDn,SDn_x,SDn_y,SDn_z,Nx,Ny,Nz);
//...........................................................................
// Gradient of the phase indicator field
//...........................................................................
Dm.CommunicateMeshHalo(SDn_x);
Dm->CommunicateMeshHalo(SDn_x);
//...........................................................................
Dm.CommunicateMeshHalo(SDn_y);
Dm->CommunicateMeshHalo(SDn_y);
//...........................................................................
Dm.CommunicateMeshHalo(SDn_z);
Dm->CommunicateMeshHalo(SDn_z);
//...........................................................................
Dm.CommunicateMeshHalo(SDs);
Dm->CommunicateMeshHalo(SDs);
pmmc_MeshGradient(SDs,SDs_x,SDs_y,SDs_z,Nx,Ny,Nz);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_x);
Dm->CommunicateMeshHalo(SDs_x);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_y);
Dm->CommunicateMeshHalo(SDs_y);
//...........................................................................
Dm.CommunicateMeshHalo(SDs_z);
Dm->CommunicateMeshHalo(SDs_z);
//...........................................................................
// Compute the mesh curvature of the phase indicator field
pmmc_MeshCurvature(SDn, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
@@ -375,26 +375,26 @@ void TwoPhase::UpdateMeshValues()
// Map Phase_tplus and Phase_tminus
for (int n=0; n<Nx*Ny*Nz; n++) dPdt(n) = 0.125*(Phase_tplus(n) - Phase_tminus(n));
//...........................................................................
Dm.CommunicateMeshHalo(Press);
Dm->CommunicateMeshHalo(Press);
//...........................................................................
Dm.CommunicateMeshHalo(Vel_x);
Dm->CommunicateMeshHalo(Vel_x);
//...........................................................................
Dm.CommunicateMeshHalo(Vel_y);
Dm->CommunicateMeshHalo(Vel_y);
//...........................................................................
Dm.CommunicateMeshHalo(Vel_z);
Dm->CommunicateMeshHalo(Vel_z);
//...........................................................................
Dm.CommunicateMeshHalo(MeanCurvature);
Dm->CommunicateMeshHalo(MeanCurvature);
//...........................................................................
Dm.CommunicateMeshHalo(GaussCurvature);
Dm->CommunicateMeshHalo(GaussCurvature);
//...........................................................................
Dm.CommunicateMeshHalo(DelPhi);
Dm->CommunicateMeshHalo(DelPhi);
//...........................................................................
// Initializing the blob ID
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;
if (Dm.id[n] == 0){
if (Dm->id[n] == 0){
// Solid phase
PhaseID(i,j,k) = 0;
}
@@ -418,8 +418,8 @@ void TwoPhase::ComputeLocal()
// If external boundary conditions are set, do not average over the inlet
kmin=1; kmax=Nz-1;
if (Dm.BoundaryCondition > 0 && Dm.kproc() == 0) kmin=4;
if (Dm.BoundaryCondition > 0 && Dm.kproc() == Dm.nprocz()-1) kmax=Nz-4;
if (Dm->BoundaryCondition > 0 && Dm->kproc() == 0) kmin=4;
if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
for (k=kmin; k<kmax; k++){
for (j=1; j<Ny-1; j++){
@@ -431,7 +431,7 @@ void TwoPhase::ComputeLocal()
// Compute volume averages
for (int p=0;p<8;p++){
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
if ( Dm.id[n] != 0 ){
if ( Dm->id[n] != 0 ){
// 1-D index for this cube corner
// compute the norm of the gradient of the phase indicator field
// Compute the non-wetting phase volume contribution
@@ -551,7 +551,7 @@ void TwoPhase::AssignComponentLabels()
int LabelNWP=1;
int LabelWP=2;
// NOTE: labeling the wetting phase components is tricky! One sandstone media had over 800,000 components
// NumberComponents_WP = ComputeGlobalPhaseComponent(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,PhaseID,LabelWP,Label_WP);
// NumberComponents_WP = ComputeGlobalPhaseComponent(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->rank_info,PhaseID,LabelWP,Label_WP);
// treat all wetting phase is connected
NumberComponents_WP=1;
for (int k=0; k<Nz; k++){
@@ -566,8 +566,8 @@ void TwoPhase::AssignComponentLabels()
}
// Fewer non-wetting phase features are present
//NumberComponents_NWP = ComputeGlobalPhaseComponent(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,PhaseID,LabelNWP,Label_NWP);
NumberComponents_NWP = ComputeGlobalBlobIDs(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,SDs,SDn,solid_isovalue,fluid_isovalue,Label_NWP,Dm.Comm);
//NumberComponents_NWP = ComputeGlobalPhaseComponent(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->rank_info,PhaseID,LabelNWP,Label_NWP);
NumberComponents_NWP = ComputeGlobalBlobIDs(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->rank_info,SDs,SDn,solid_isovalue,fluid_isovalue,Label_NWP,Dm->Comm);
}
void TwoPhase::ComponentAverages()
@@ -585,15 +585,15 @@ void TwoPhase::ComponentAverages()
ComponentAverages_WP.fill(0.0);
ComponentAverages_NWP.fill(0.0);
if (Dm.rank()==0){
if (Dm->rank()==0){
printf("Number of wetting phase components is %i \n",NumberComponents_WP);
printf("Number of non-wetting phase components is %i \n",NumberComponents_NWP);
}
// If external boundary conditions are set, do not average over the inlet
kmin=1; kmax=Nz-1;
if (Dm.BoundaryCondition > 0 && Dm.kproc() == 0) kmin=4;
if (Dm.BoundaryCondition > 0 && Dm.kproc() == Dm.nprocz()-1) kmax=Nz-4;
if (Dm->BoundaryCondition > 0 && Dm->kproc() == 0) kmin=4;
if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
for (k=kmin; k<kmax; k++){
for (j=1; j<Ny-1; j++){
@@ -624,7 +624,7 @@ void TwoPhase::ComponentAverages()
// Compute volume averages
for (int p=0;p<8;p++){
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
if ( Dm.id[n] != 0 ){
if ( Dm->id[n] != 0 ){
// 1-D index for this cube corner
// compute the norm of the gradient of the phase indicator field
// Compute the non-wetting phase volume contribution
@@ -636,9 +636,9 @@ void TwoPhase::ComponentAverages()
ComponentAverages_NWP(VY,LabelNWP) += 0.125*Vel_y(n);
ComponentAverages_NWP(VZ,LabelNWP) += 0.125*Vel_z(n);
// center of mass
ComponentAverages_NWP(CMX,LabelNWP) += 0.125*(i+cube[p][0]+Dm.iproc()*Nx);
ComponentAverages_NWP(CMY,LabelNWP) += 0.125*(j+cube[p][1]+Dm.jproc()*Ny);
ComponentAverages_NWP(CMZ,LabelNWP) += 0.125*(k+cube[p][2]+Dm.kproc()*Nz);
ComponentAverages_NWP(CMX,LabelNWP) += 0.125*(i+cube[p][0]+Dm->iproc()*Nx);
ComponentAverages_NWP(CMY,LabelNWP) += 0.125*(j+cube[p][1]+Dm->jproc()*Ny);
ComponentAverages_NWP(CMZ,LabelNWP) += 0.125*(k+cube[p][2]+Dm->kproc()*Nz);
// twice the kinetic energy
ComponentAverages_NWP(VSQ,LabelNWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n));
@@ -656,9 +656,9 @@ void TwoPhase::ComponentAverages()
ComponentAverages_WP(VY,LabelWP)+= 0.125*Vel_y(n);
ComponentAverages_WP(VZ,LabelWP) += 0.125*Vel_z(n);
// Center of mass
ComponentAverages_WP(CMX,LabelWP) += 0.125*(i+cube[p][0]+Dm.iproc()*Nx);
ComponentAverages_WP(CMY,LabelWP) += 0.125*(j+cube[p][1]+Dm.jproc()*Ny);
ComponentAverages_WP(CMZ,LabelWP) += 0.125*(k+cube[p][2]+Dm.kproc()*Nz);
ComponentAverages_WP(CMX,LabelWP) += 0.125*(i+cube[p][0]+Dm->iproc()*Nx);
ComponentAverages_WP(CMY,LabelWP) += 0.125*(j+cube[p][1]+Dm->jproc()*Ny);
ComponentAverages_WP(CMZ,LabelWP) += 0.125*(k+cube[p][2]+Dm->kproc()*Nz);
// twice the kinetic energy
ComponentAverages_WP(VSQ,LabelWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n));
@@ -802,24 +802,24 @@ void TwoPhase::ComponentAverages()
}
}
MPI_Barrier(Dm.Comm);
if (Dm.rank()==0){
MPI_Barrier(Dm->Comm);
if (Dm->rank()==0){
printf("Component averages computed locally -- reducing result... \n");
}
// Globally reduce the non-wetting phase averages
RecvBuffer.resize(BLOB_AVG_COUNT,NumberComponents_NWP);
/* for (int b=0; b<NumberComponents_NWP; b++){
MPI_Barrier(Dm.Comm);
MPI_Allreduce(&ComponentAverages_NWP(0,b),&RecvBuffer(0),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Barrier(Dm->Comm);
MPI_Allreduce(&ComponentAverages_NWP(0,b),&RecvBuffer(0),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,Dm->Comm);
for (int idx=0; idx<BLOB_AVG_COUNT; idx++) ComponentAverages_NWP(idx,b)=RecvBuffer(idx);
}
*/
MPI_Barrier(Dm.Comm);
MPI_Allreduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT*NumberComponents_NWP, MPI_DOUBLE,MPI_SUM,Dm.Comm);
// MPI_Reduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
MPI_Barrier(Dm->Comm);
MPI_Allreduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT*NumberComponents_NWP, MPI_DOUBLE,MPI_SUM,Dm->Comm);
// MPI_Reduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm->Comm);
if (Dm.rank()==0){
if (Dm->rank()==0){
printf("rescaling... \n");
}
@@ -907,15 +907,15 @@ void TwoPhase::ComponentAverages()
}
}
if (Dm.rank()==0){
if (Dm->rank()==0){
printf("reduce WP averages... \n");
}
// reduce the wetting phase averages
for (int b=0; b<NumberComponents_WP; b++){
MPI_Barrier(Dm.Comm);
// MPI_Allreduce(&ComponentAverages_WP(0,b),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Reduce(&ComponentAverages_WP(0,b),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
MPI_Barrier(Dm->Comm);
// MPI_Allreduce(&ComponentAverages_WP(0,b),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Reduce(&ComponentAverages_WP(0,b),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm->Comm);
for (int idx=0; idx<BLOB_AVG_COUNT; idx++) ComponentAverages_WP(idx,b)=RecvBuffer(idx);
}
@@ -1052,15 +1052,15 @@ void TwoPhase::WriteSurfaces(int logcount)
C = P;
}
// Remap the points
A.x += 1.0*Dm.iproc()*(Nx-2);
A.y += 1.0*Dm.jproc()*(Nx-2);
A.z += 1.0*Dm.kproc()*(Nx-2);
B.x += 1.0*Dm.iproc()*(Nx-2);
B.y += 1.0*Dm.jproc()*(Nx-2);
B.z += 1.0*Dm.kproc()*(Nx-2);
C.x += 1.0*Dm.iproc()*(Nx-2);
C.y += 1.0*Dm.jproc()*(Nx-2);
C.z += 1.0*Dm.kproc()*(Nx-2);
A.x += 1.0*Dm->iproc()*(Nx-2);
A.y += 1.0*Dm->jproc()*(Nx-2);
A.z += 1.0*Dm->kproc()*(Nx-2);
B.x += 1.0*Dm->iproc()*(Nx-2);
B.y += 1.0*Dm->jproc()*(Nx-2);
B.z += 1.0*Dm->kproc()*(Nx-2);
C.x += 1.0*Dm->iproc()*(Nx-2);
C.y += 1.0*Dm->jproc()*(Nx-2);
C.z += 1.0*Dm->kproc()*(Nx-2);
wn_mesh->A.push_back(A);
wn_mesh->B.push_back(B);
wn_mesh->C.push_back(C);
@@ -1070,15 +1070,15 @@ void TwoPhase::WriteSurfaces(int logcount)
B = ws_pts(ws_tris(1,r));
C = ws_pts(ws_tris(2,r));
// Remap the points
A.x += 1.0*Dm.iproc()*(Nx-2);
A.y += 1.0*Dm.jproc()*(Nx-2);
A.z += 1.0*Dm.kproc()*(Nx-2);
B.x += 1.0*Dm.iproc()*(Nx-2);
B.y += 1.0*Dm.jproc()*(Nx-2);
B.z += 1.0*Dm.kproc()*(Nx-2);
C.x += 1.0*Dm.iproc()*(Nx-2);
C.y += 1.0*Dm.jproc()*(Nx-2);
C.z += 1.0*Dm.kproc()*(Nx-2);
A.x += 1.0*Dm->iproc()*(Nx-2);
A.y += 1.0*Dm->jproc()*(Nx-2);
A.z += 1.0*Dm->kproc()*(Nx-2);
B.x += 1.0*Dm->iproc()*(Nx-2);
B.y += 1.0*Dm->jproc()*(Nx-2);
B.z += 1.0*Dm->kproc()*(Nx-2);
C.x += 1.0*Dm->iproc()*(Nx-2);
C.y += 1.0*Dm->jproc()*(Nx-2);
C.z += 1.0*Dm->kproc()*(Nx-2);
ws_mesh->A.push_back(A);
ws_mesh->B.push_back(B);
ws_mesh->C.push_back(C);
@@ -1088,15 +1088,15 @@ void TwoPhase::WriteSurfaces(int logcount)
B = ns_pts(ns_tris(1,r));
C = ns_pts(ns_tris(2,r));
// Remap the points
A.x += 1.0*Dm.iproc()*(Nx-2);
A.y += 1.0*Dm.jproc()*(Nx-2);
A.z += 1.0*Dm.kproc()*(Nx-2);
B.x += 1.0*Dm.iproc()*(Nx-2);
B.y += 1.0*Dm.jproc()*(Nx-2);
B.z += 1.0*Dm.kproc()*(Nx-2);
C.x += 1.0*Dm.iproc()*(Nx-2);
C.y += 1.0*Dm.jproc()*(Nx-2);
C.z += 1.0*Dm.kproc()*(Nx-2);
A.x += 1.0*Dm->iproc()*(Nx-2);
A.y += 1.0*Dm->jproc()*(Nx-2);
A.z += 1.0*Dm->kproc()*(Nx-2);
B.x += 1.0*Dm->iproc()*(Nx-2);
B.y += 1.0*Dm->jproc()*(Nx-2);
B.z += 1.0*Dm->kproc()*(Nx-2);
C.x += 1.0*Dm->iproc()*(Nx-2);
C.y += 1.0*Dm->jproc()*(Nx-2);
C.z += 1.0*Dm->kproc()*(Nx-2);
ns_mesh->A.push_back(A);
ns_mesh->B.push_back(B);
ns_mesh->C.push_back(C);
@@ -1112,7 +1112,7 @@ void TwoPhase::WriteSurfaces(int logcount)
meshData[1].mesh = ws_mesh;
meshData[2].meshName = "ns-tris";
meshData[2].mesh = ns_mesh;
IO::writeData( logcount, meshData, Dm.Comm );
IO::writeData( logcount, meshData, Dm->Comm );
}
@@ -1121,43 +1121,43 @@ void TwoPhase::Reduce()
int i;
double iVol_global=1.0/Volume;
//...........................................................................
MPI_Barrier(Dm.Comm);
MPI_Allreduce(&nwp_volume,&nwp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&wp_volume,&wp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&awn,&awn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&ans,&ans_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&aws,&aws_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&lwns,&lwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&As,&As_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Jwn,&Jwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Kwn,&Kwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&KGwns,&KGwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&KNwns,&KNwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&wwndnw,&wwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&wwnsdnwn,&wwnsdnwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Jwnwwndnw,&Jwnwwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Barrier(Dm->Comm);
MPI_Allreduce(&nwp_volume,&nwp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&wp_volume,&wp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&awn,&awn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&ans,&ans_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&aws,&aws_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&lwns,&lwns_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&As,&As_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&Jwn,&Jwn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&Kwn,&Kwn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&KGwns,&KGwns_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&KNwns,&KNwns_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&wwndnw,&wwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&wwnsdnwn,&wwnsdnwn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&Jwnwwndnw,&Jwnwwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
// Phase averages
MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&paw,&paw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&pan,&pan_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vaw(0),&vaw_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&van(0),&van_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vawn(0),&vawn_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vawns(0),&vawns_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Gwn(0),&Gwn_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Gns(0),&Gns_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Gws(0),&Gws_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&trawn,&trawn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&trJwn,&trJwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&trRwn,&trRwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&euler,&euler_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&An,&An_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Jn,&Jn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Kn,&Kn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&paw,&paw_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&pan,&pan_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&vaw(0),&vaw_global(0),3,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&van(0),&van_global(0),3,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&vawn(0),&vawn_global(0),3,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&vawns(0),&vawns_global(0),3,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&Gwn(0),&Gwn_global(0),6,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&Gns(0),&Gns_global(0),6,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&Gws(0),&Gws_global(0),6,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&trawn,&trawn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&trJwn,&trJwn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&trRwn,&trRwn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&euler,&euler_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&An,&An_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&Jn,&Jn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&Kn,&Kn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Barrier(Dm.Comm);
MPI_Barrier(Dm->Comm);
// Normalize the phase averages
// (density of both components = 1.0)
@@ -1224,7 +1224,7 @@ void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT)
void TwoPhase::PrintAll(int timestep)
{
if (Dm.rank()==0){
if (Dm->rank()==0){
fprintf(TIMELOG,"%i %.5g ",timestep,dEs); // change in surface energy
fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure
fprintf(TIMELOG,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas
@@ -1277,7 +1277,7 @@ void TwoPhase::PrintAll(int timestep)
void TwoPhase::PrintComponents(int timestep)
{
if (Dm.rank()==0){
if (Dm->rank()==0){
printf("PRINT %i COMPONENT AVEREAGES: time = %i \n",(int)ComponentAverages_NWP.size(1),timestep);
for (int b=0; b<NumberComponents_NWP; b++){
//if (ComponentAverages_NWP(TRIMVOL,b) > 0.0){

4
analysis/TwoPhase.h
View File

@@ -64,7 +64,7 @@ class TwoPhase{
public:
//...........................................................................
Domain& Dm;
Domain* Dm;
int NumberComponents_WP,NumberComponents_NWP;
//...........................................................................
// Averaging variables
@@ -145,7 +145,7 @@ public:
DoubleArray ComponentAverages_WP;
DoubleArray ComponentAverages_NWP;
//...........................................................................
TwoPhase(Domain &dm);
TwoPhase(Domain *dm);
~TwoPhase();
void Initialize();
// void SetupCubes(Domain &Dm);

24
analysis/runAnalysis.cpp
View File

@@ -282,22 +282,22 @@ runAnalysis::commWrapper runAnalysis::getComm( )
* Constructor/Destructors *
******************************************************************/
runAnalysis::runAnalysis( std::shared_ptr<Database> db,
const RankInfoStruct& rank_info, const ScaLBL_Communicator &ScaLBL_Comm, const Domain& Dm,
const RankInfoStruct& rank_info, const ScaLBL_Communicator &ScaLBL_Comm, const Domain* Dm,
int Np, bool pBC, double beta, IntArray Map ):
d_Np( Np ),
d_beta( beta ),
d_ScaLBL_Comm( ScaLBL_Comm ),
d_rank_info( rank_info ),
d_Map( Map ),
d_fillData(Dm.Comm,Dm.rank_info,Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,1,1,1,0,1)
d_ScaLBL_Comm( ScaLBL_Comm ),
d_fillData(Dm->Comm,Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,1,1,1,0,1)
{
NULL_USE( pBC );
INSIST( db, "Input database is empty" );
char rankString[20];
sprintf(rankString,"%05d",Dm.rank());
d_N[0] = Dm.Nx;
d_N[1] = Dm.Ny;
d_N[2] = Dm.Nz;
sprintf(rankString,"%05d",Dm->rank());
d_N[0] = Dm->Nx;
d_N[1] = Dm->Ny;
d_N[2] = Dm->Nz;
d_restart_interval = db->getScalar<int>( "restart_interval" );
d_analysis_interval = db->getScalar<int>( "analysis_interval" );
d_blobid_interval = db->getScalar<int>( "blobid_interval" );
@@ -309,7 +309,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
// Create the MeshDataStruct
d_meshData.resize(1);
d_meshData[0].meshName = "domain";
d_meshData[0].mesh = std::make_shared<IO::DomainMesh>( Dm.rank_info,Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.Lx,Dm.Ly,Dm.Lz );
d_meshData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
auto PhaseVar = std::make_shared<IO::Variable>();
auto PressVar = std::make_shared<IO::Variable>();
auto SignDistVar = std::make_shared<IO::Variable>();
@@ -317,22 +317,22 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
PhaseVar->name = "phase";
PhaseVar->type = IO::VariableType::VolumeVariable;
PhaseVar->dim = 1;
PhaseVar->data.resize(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2);
PhaseVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
d_meshData[0].vars.push_back(PhaseVar);
PressVar->name = "Pressure";
PressVar->type = IO::VariableType::VolumeVariable;
PressVar->dim = 1;
PressVar->data.resize(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2);
PressVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
d_meshData[0].vars.push_back(PressVar);
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VariableType::VolumeVariable;
SignDistVar->dim = 1;
SignDistVar->data.resize(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2);
SignDistVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
d_meshData[0].vars.push_back(SignDistVar);
BlobIDVar->name = "BlobID";
BlobIDVar->type = IO::VariableType::VolumeVariable;
BlobIDVar->dim = 1;
BlobIDVar->data.resize(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2);
BlobIDVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
d_meshData[0].vars.push_back(BlobIDVar);
// Initialize the comms
MPI_Comm_dup(MPI_COMM_WORLD,&d_comm);

4
analysis/runAnalysis.h
View File

@@ -24,13 +24,13 @@ public:
//! Constructor
runAnalysis( std::shared_ptr<Database> db, const RankInfoStruct& rank_info,
const ScaLBL_Communicator &ScaLBL_Comm, const Domain& dm, int Np, bool pBC, double beta, IntArray Map );
const ScaLBL_Communicator &ScaLBL_Comm, const Domain* dm, int Np, bool pBC, double beta, IntArray Map );
//! Destructor
~runAnalysis();
//! Run the next analysis
void run( int timestep, TwoPhase& Averages, const double *Phi,
void run( int timestep, TwoPhase &Averages, const double *Phi,
double *Pressure, double *Velocity, double *fq, double *Den );
//! Finish all active analysis

294
common/ScaLBL.cpp
View File

@@ -1,83 +1,83 @@
#include "common/ScaLBL.h"
ScaLBL_Communicator::ScaLBL_Communicator(Domain &Dm){
ScaLBL_Communicator::ScaLBL_Communicator(Domain *Dm){
//......................................................................................
Lock=false; // unlock the communicator
//......................................................................................
// Create a separate copy of the communicator for the device
//MPI_Comm_group(Dm.Comm,&Group);
//MPI_Comm_create(Dm.Comm,Group,&MPI_COMM_SCALBL);
MPI_Comm_dup(Dm.Comm,&MPI_COMM_SCALBL);
//MPI_Comm_group(Dm->Comm,&Group);
//MPI_Comm_create(Dm->Comm,Group,&MPI_COMM_SCALBL);
MPI_Comm_dup(Dm->Comm,&MPI_COMM_SCALBL);
//......................................................................................
// Copy the domain size and communication information directly from Dm
Nx = Dm.Nx;
Ny = Dm.Ny;
Nz = Dm.Nz;
Nx = Dm->Nx;
Ny = Dm->Ny;
Nz = Dm->Nz;
N = Nx*Ny*Nz;
next=0;
rank=Dm.rank();
rank_x=Dm.rank_x();
rank_y=Dm.rank_y();
rank_z=Dm.rank_z();
rank_X=Dm.rank_X();
rank_Y=Dm.rank_Y();
rank_Z=Dm.rank_Z();
rank_xy=Dm.rank_xy();
rank_XY=Dm.rank_XY();
rank_xY=Dm.rank_xY();
rank_Xy=Dm.rank_Xy();
rank_xz=Dm.rank_xz();
rank_XZ=Dm.rank_XZ();
rank_xZ=Dm.rank_xZ();
rank_Xz=Dm.rank_Xz();
rank_yz=Dm.rank_yz();
rank_YZ=Dm.rank_YZ();
rank_yZ=Dm.rank_yZ();
rank_Yz=Dm.rank_Yz();
sendCount_x=Dm.sendCount_x;
sendCount_y=Dm.sendCount_y;
sendCount_z=Dm.sendCount_z;
sendCount_X=Dm.sendCount_X;
sendCount_Y=Dm.sendCount_Y;
sendCount_Z=Dm.sendCount_Z;
sendCount_xy=Dm.sendCount_xy;
sendCount_yz=Dm.sendCount_yz;
sendCount_xz=Dm.sendCount_xz;
sendCount_Xy=Dm.sendCount_Xy;
sendCount_Yz=Dm.sendCount_Yz;
sendCount_xZ=Dm.sendCount_xZ;
sendCount_xY=Dm.sendCount_xY;
sendCount_yZ=Dm.sendCount_yZ;
sendCount_Xz=Dm.sendCount_Xz;
sendCount_XY=Dm.sendCount_XY;
sendCount_YZ=Dm.sendCount_YZ;
sendCount_XZ=Dm.sendCount_XZ;
recvCount_x=Dm.recvCount_x;
recvCount_y=Dm.recvCount_y;
recvCount_z=Dm.recvCount_z;
recvCount_X=Dm.recvCount_X;
recvCount_Y=Dm.recvCount_Y;
recvCount_Z=Dm.recvCount_Z;
recvCount_xy=Dm.recvCount_xy;
recvCount_yz=Dm.recvCount_yz;
recvCount_xz=Dm.recvCount_xz;
recvCount_Xy=Dm.recvCount_Xy;
recvCount_Yz=Dm.recvCount_Yz;
recvCount_xZ=Dm.recvCount_xZ;
recvCount_xY=Dm.recvCount_xY;
recvCount_yZ=Dm.recvCount_yZ;
recvCount_Xz=Dm.recvCount_Xz;
recvCount_XY=Dm.recvCount_XY;
recvCount_YZ=Dm.recvCount_YZ;
recvCount_XZ=Dm.recvCount_XZ;
rank=Dm->rank();
rank_x=Dm->rank_x();
rank_y=Dm->rank_y();
rank_z=Dm->rank_z();
rank_X=Dm->rank_X();
rank_Y=Dm->rank_Y();
rank_Z=Dm->rank_Z();
rank_xy=Dm->rank_xy();
rank_XY=Dm->rank_XY();
rank_xY=Dm->rank_xY();
rank_Xy=Dm->rank_Xy();
rank_xz=Dm->rank_xz();
rank_XZ=Dm->rank_XZ();
rank_xZ=Dm->rank_xZ();
rank_Xz=Dm->rank_Xz();
rank_yz=Dm->rank_yz();
rank_YZ=Dm->rank_YZ();
rank_yZ=Dm->rank_yZ();
rank_Yz=Dm->rank_Yz();
sendCount_x=Dm->sendCount_x;
sendCount_y=Dm->sendCount_y;
sendCount_z=Dm->sendCount_z;
sendCount_X=Dm->sendCount_X;
sendCount_Y=Dm->sendCount_Y;
sendCount_Z=Dm->sendCount_Z;
sendCount_xy=Dm->sendCount_xy;
sendCount_yz=Dm->sendCount_yz;
sendCount_xz=Dm->sendCount_xz;
sendCount_Xy=Dm->sendCount_Xy;
sendCount_Yz=Dm->sendCount_Yz;
sendCount_xZ=Dm->sendCount_xZ;
sendCount_xY=Dm->sendCount_xY;
sendCount_yZ=Dm->sendCount_yZ;
sendCount_Xz=Dm->sendCount_Xz;
sendCount_XY=Dm->sendCount_XY;
sendCount_YZ=Dm->sendCount_YZ;
sendCount_XZ=Dm->sendCount_XZ;
recvCount_x=Dm->recvCount_x;
recvCount_y=Dm->recvCount_y;
recvCount_z=Dm->recvCount_z;
recvCount_X=Dm->recvCount_X;
recvCount_Y=Dm->recvCount_Y;
recvCount_Z=Dm->recvCount_Z;
recvCount_xy=Dm->recvCount_xy;
recvCount_yz=Dm->recvCount_yz;
recvCount_xz=Dm->recvCount_xz;
recvCount_Xy=Dm->recvCount_Xy;
recvCount_Yz=Dm->recvCount_Yz;
recvCount_xZ=Dm->recvCount_xZ;
recvCount_xY=Dm->recvCount_xY;
recvCount_yZ=Dm->recvCount_yZ;
recvCount_Xz=Dm->recvCount_Xz;
recvCount_XY=Dm->recvCount_XY;
recvCount_YZ=Dm->recvCount_YZ;
recvCount_XZ=Dm->recvCount_XZ;
iproc = Dm.iproc();
jproc = Dm.jproc();
kproc = Dm.kproc();
nprocx = Dm.nprocx();
nprocy = Dm.nprocy();
nprocz = Dm.nprocz();
BoundaryCondition = Dm.BoundaryCondition;
iproc = Dm->iproc();
jproc = Dm->jproc();
kproc = Dm->kproc();
nprocx = Dm->nprocx();
nprocy = Dm->nprocy();
nprocz = Dm->nprocz();
BoundaryCondition = Dm->BoundaryCondition;
//......................................................................................
ScaLBL_AllocateZeroCopy((void **) &sendbuf_x, 5*sendCount_x*sizeof(double)); // Allocate device memory
@@ -176,43 +176,43 @@ ScaLBL_Communicator::ScaLBL_Communicator(Domain &Dm){
ScaLBL_AllocateZeroCopy((void **) &dvcRecvDist_YZ, recvCount_YZ*sizeof(int)); // Allocate device memory
//......................................................................................
ScaLBL_CopyToZeroCopy(dvcSendList_x,Dm.sendList_x,sendCount_x*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_X,Dm.sendList_X,sendCount_X*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_y,Dm.sendList_y,sendCount_y*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Y,Dm.sendList_Y,sendCount_Y*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_z,Dm.sendList_z,sendCount_z*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Z,Dm.sendList_Z,sendCount_Z*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_xy,Dm.sendList_xy,sendCount_xy*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_XY,Dm.sendList_XY,sendCount_XY*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_xY,Dm.sendList_xY,sendCount_xY*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Xy,Dm.sendList_Xy,sendCount_Xy*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_xz,Dm.sendList_xz,sendCount_xz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_XZ,Dm.sendList_XZ,sendCount_XZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_xZ,Dm.sendList_xZ,sendCount_xZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Xz,Dm.sendList_Xz,sendCount_Xz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_yz,Dm.sendList_yz,sendCount_yz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_YZ,Dm.sendList_YZ,sendCount_YZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_yZ,Dm.sendList_yZ,sendCount_yZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Yz,Dm.sendList_Yz,sendCount_Yz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_x,Dm->sendList_x,sendCount_x*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_X,Dm->sendList_X,sendCount_X*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_y,Dm->sendList_y,sendCount_y*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Y,Dm->sendList_Y,sendCount_Y*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_z,Dm->sendList_z,sendCount_z*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Z,Dm->sendList_Z,sendCount_Z*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_xy,Dm->sendList_xy,sendCount_xy*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_XY,Dm->sendList_XY,sendCount_XY*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_xY,Dm->sendList_xY,sendCount_xY*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Xy,Dm->sendList_Xy,sendCount_Xy*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_xz,Dm->sendList_xz,sendCount_xz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_XZ,Dm->sendList_XZ,sendCount_XZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_xZ,Dm->sendList_xZ,sendCount_xZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Xz,Dm->sendList_Xz,sendCount_Xz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_yz,Dm->sendList_yz,sendCount_yz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_YZ,Dm->sendList_YZ,sendCount_YZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_yZ,Dm->sendList_yZ,sendCount_yZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcSendList_Yz,Dm->sendList_Yz,sendCount_Yz*sizeof(int));
//......................................................................................
ScaLBL_CopyToZeroCopy(dvcRecvList_x,Dm.recvList_x,recvCount_x*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_X,Dm.recvList_X,recvCount_X*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_y,Dm.recvList_y,recvCount_y*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Y,Dm.recvList_Y,recvCount_Y*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_z,Dm.recvList_z,recvCount_z*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Z,Dm.recvList_Z,recvCount_Z*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_xy,Dm.recvList_xy,recvCount_xy*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_XY,Dm.recvList_XY,recvCount_XY*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_xY,Dm.recvList_xY,recvCount_xY*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Xy,Dm.recvList_Xy,recvCount_Xy*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_xz,Dm.recvList_xz,recvCount_xz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_XZ,Dm.recvList_XZ,recvCount_XZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_xZ,Dm.recvList_xZ,recvCount_xZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Xz,Dm.recvList_Xz,recvCount_Xz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_yz,Dm.recvList_yz,recvCount_yz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_YZ,Dm.recvList_YZ,recvCount_YZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_yZ,Dm.recvList_yZ,recvCount_yZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Yz,Dm.recvList_Yz,recvCount_Yz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_x,Dm->recvList_x,recvCount_x*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_X,Dm->recvList_X,recvCount_X*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_y,Dm->recvList_y,recvCount_y*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Y,Dm->recvList_Y,recvCount_Y*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_z,Dm->recvList_z,recvCount_z*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Z,Dm->recvList_Z,recvCount_Z*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_xy,Dm->recvList_xy,recvCount_xy*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_XY,Dm->recvList_XY,recvCount_XY*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_xY,Dm->recvList_xY,recvCount_xY*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Xy,Dm->recvList_Xy,recvCount_Xy*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_xz,Dm->recvList_xz,recvCount_xz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_XZ,Dm->recvList_XZ,recvCount_XZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_xZ,Dm->recvList_xZ,recvCount_xZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Xz,Dm->recvList_Xz,recvCount_Xz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_yz,Dm->recvList_yz,recvCount_yz*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_YZ,Dm->recvList_YZ,recvCount_YZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_yZ,Dm->recvList_yZ,recvCount_yZ*sizeof(int));
ScaLBL_CopyToZeroCopy(dvcRecvList_Yz,Dm->recvList_Yz,recvCount_Yz*sizeof(int));
//......................................................................................
MPI_Barrier(MPI_COMM_SCALBL);
@@ -221,70 +221,70 @@ ScaLBL_Communicator::ScaLBL_Communicator(Domain &Dm){
// Set up the recieve distribution lists
//...................................................................................
//...Map recieve list for the X face: q=2,8,10,12,14 .................................
D3Q19_MapRecv(-1,0,0,Dm.recvList_X,0,recvCount_X,dvcRecvDist_X);
D3Q19_MapRecv(-1,-1,0,Dm.recvList_X,recvCount_X,recvCount_X,dvcRecvDist_X);
D3Q19_MapRecv(-1,1,0,Dm.recvList_X,2*recvCount_X,recvCount_X,dvcRecvDist_X);
D3Q19_MapRecv(-1,0,-1,Dm.recvList_X,3*recvCount_X,recvCount_X,dvcRecvDist_X);
D3Q19_MapRecv(-1,0,1,Dm.recvList_X,4*recvCount_X,recvCount_X,dvcRecvDist_X);
D3Q19_MapRecv(-1,0,0,Dm->recvList_X,0,recvCount_X,dvcRecvDist_X);
D3Q19_MapRecv(-1,-1,0,Dm->recvList_X,recvCount_X,recvCount_X,dvcRecvDist_X);
D3Q19_MapRecv(-1,1,0,Dm->recvList_X,2*recvCount_X,recvCount_X,dvcRecvDist_X);
D3Q19_MapRecv(-1,0,-1,Dm->recvList_X,3*recvCount_X,recvCount_X,dvcRecvDist_X);
D3Q19_MapRecv(-1,0,1,Dm->recvList_X,4*recvCount_X,recvCount_X,dvcRecvDist_X);
//...................................................................................
//...Map recieve list for the x face: q=1,7,9,11,13..................................
D3Q19_MapRecv(1,0,0,Dm.recvList_x,0,recvCount_x,dvcRecvDist_x);
D3Q19_MapRecv(1,1,0,Dm.recvList_x,recvCount_x,recvCount_x,dvcRecvDist_x);
D3Q19_MapRecv(1,-1,0,Dm.recvList_x,2*recvCount_x,recvCount_x,dvcRecvDist_x);
D3Q19_MapRecv(1,0,1,Dm.recvList_x,3*recvCount_x,recvCount_x,dvcRecvDist_x);
D3Q19_MapRecv(1,0,-1,Dm.recvList_x,4*recvCount_x,recvCount_x,dvcRecvDist_x);
D3Q19_MapRecv(1,0,0,Dm->recvList_x,0,recvCount_x,dvcRecvDist_x);
D3Q19_MapRecv(1,1,0,Dm->recvList_x,recvCount_x,recvCount_x,dvcRecvDist_x);
D3Q19_MapRecv(1,-1,0,Dm->recvList_x,2*recvCount_x,recvCount_x,dvcRecvDist_x);
D3Q19_MapRecv(1,0,1,Dm->recvList_x,3*recvCount_x,recvCount_x,dvcRecvDist_x);
D3Q19_MapRecv(1,0,-1,Dm->recvList_x,4*recvCount_x,recvCount_x,dvcRecvDist_x);
//...................................................................................
//...Map recieve list for the y face: q=4,8,9,16,18 ...................................
D3Q19_MapRecv(0,-1,0,Dm.recvList_Y,0,recvCount_Y,dvcRecvDist_Y);
D3Q19_MapRecv(-1,-1,0,Dm.recvList_Y,recvCount_Y,recvCount_Y,dvcRecvDist_Y);
D3Q19_MapRecv(1,-1,0,Dm.recvList_Y,2*recvCount_Y,recvCount_Y,dvcRecvDist_Y);
D3Q19_MapRecv(0,-1,-1,Dm.recvList_Y,3*recvCount_Y,recvCount_Y,dvcRecvDist_Y);
D3Q19_MapRecv(0,-1,1,Dm.recvList_Y,4*recvCount_Y,recvCount_Y,dvcRecvDist_Y);
D3Q19_MapRecv(0,-1,0,Dm->recvList_Y,0,recvCount_Y,dvcRecvDist_Y);
D3Q19_MapRecv(-1,-1,0,Dm->recvList_Y,recvCount_Y,recvCount_Y,dvcRecvDist_Y);
D3Q19_MapRecv(1,-1,0,Dm->recvList_Y,2*recvCount_Y,recvCount_Y,dvcRecvDist_Y);
D3Q19_MapRecv(0,-1,-1,Dm->recvList_Y,3*recvCount_Y,recvCount_Y,dvcRecvDist_Y);
D3Q19_MapRecv(0,-1,1,Dm->recvList_Y,4*recvCount_Y,recvCount_Y,dvcRecvDist_Y);
//...................................................................................
//...Map recieve list for the Y face: q=3,7,10,15,17 ..................................
D3Q19_MapRecv(0,1,0,Dm.recvList_y,0,recvCount_y,dvcRecvDist_y);
D3Q19_MapRecv(1,1,0,Dm.recvList_y,recvCount_y,recvCount_y,dvcRecvDist_y);
D3Q19_MapRecv(-1,1,0,Dm.recvList_y,2*recvCount_y,recvCount_y,dvcRecvDist_y);
D3Q19_MapRecv(0,1,1,Dm.recvList_y,3*recvCount_y,recvCount_y,dvcRecvDist_y);
D3Q19_MapRecv(0,1,-1,Dm.recvList_y,4*recvCount_y,recvCount_y,dvcRecvDist_y);
D3Q19_MapRecv(0,1,0,Dm->recvList_y,0,recvCount_y,dvcRecvDist_y);
D3Q19_MapRecv(1,1,0,Dm->recvList_y,recvCount_y,recvCount_y,dvcRecvDist_y);
D3Q19_MapRecv(-1,1,0,Dm->recvList_y,2*recvCount_y,recvCount_y,dvcRecvDist_y);
D3Q19_MapRecv(0,1,1,Dm->recvList_y,3*recvCount_y,recvCount_y,dvcRecvDist_y);
D3Q19_MapRecv(0,1,-1,Dm->recvList_y,4*recvCount_y,recvCount_y,dvcRecvDist_y);
//...................................................................................
//...Map recieve list for the z face<<<6,12,13,16,17)..............................................
D3Q19_MapRecv(0,0,-1,Dm.recvList_Z,0,recvCount_Z,dvcRecvDist_Z);
D3Q19_MapRecv(-1,0,-1,Dm.recvList_Z,recvCount_Z,recvCount_Z,dvcRecvDist_Z);
D3Q19_MapRecv(1,0,-1,Dm.recvList_Z,2*recvCount_Z,recvCount_Z,dvcRecvDist_Z);
D3Q19_MapRecv(0,-1,-1,Dm.recvList_Z,3*recvCount_Z,recvCount_Z,dvcRecvDist_Z);
D3Q19_MapRecv(0,1,-1,Dm.recvList_Z,4*recvCount_Z,recvCount_Z,dvcRecvDist_Z);
D3Q19_MapRecv(0,0,-1,Dm->recvList_Z,0,recvCount_Z,dvcRecvDist_Z);
D3Q19_MapRecv(-1,0,-1,Dm->recvList_Z,recvCount_Z,recvCount_Z,dvcRecvDist_Z);
D3Q19_MapRecv(1,0,-1,Dm->recvList_Z,2*recvCount_Z,recvCount_Z,dvcRecvDist_Z);
D3Q19_MapRecv(0,-1,-1,Dm->recvList_Z,3*recvCount_Z,recvCount_Z,dvcRecvDist_Z);
D3Q19_MapRecv(0,1,-1,Dm->recvList_Z,4*recvCount_Z,recvCount_Z,dvcRecvDist_Z);
//...Map recieve list for the Z face<<<5,11,14,15,18)..............................................
D3Q19_MapRecv(0,0,1,Dm.recvList_z,0,recvCount_z,dvcRecvDist_z);
D3Q19_MapRecv(1,0,1,Dm.recvList_z,recvCount_z,recvCount_z,dvcRecvDist_z);
D3Q19_MapRecv(-1,0,1,Dm.recvList_z,2*recvCount_z,recvCount_z,dvcRecvDist_z);
D3Q19_MapRecv(0,1,1,Dm.recvList_z,3*recvCount_z,recvCount_z,dvcRecvDist_z);
D3Q19_MapRecv(0,-1,1,Dm.recvList_z,4*recvCount_z,recvCount_z,dvcRecvDist_z);
D3Q19_MapRecv(0,0,1,Dm->recvList_z,0,recvCount_z,dvcRecvDist_z);
D3Q19_MapRecv(1,0,1,Dm->recvList_z,recvCount_z,recvCount_z,dvcRecvDist_z);
D3Q19_MapRecv(-1,0,1,Dm->recvList_z,2*recvCount_z,recvCount_z,dvcRecvDist_z);
D3Q19_MapRecv(0,1,1,Dm->recvList_z,3*recvCount_z,recvCount_z,dvcRecvDist_z);
D3Q19_MapRecv(0,-1,1,Dm->recvList_z,4*recvCount_z,recvCount_z,dvcRecvDist_z);
//..................................................................................
//...Map recieve list for the xy edge <<<8)................................
D3Q19_MapRecv(-1,-1,0,Dm.recvList_XY,0,recvCount_XY,dvcRecvDist_XY);
D3Q19_MapRecv(-1,-1,0,Dm->recvList_XY,0,recvCount_XY,dvcRecvDist_XY);
//...Map recieve list for the Xy edge <<<9)................................
D3Q19_MapRecv(1,-1,0,Dm.recvList_xY,0,recvCount_xY,dvcRecvDist_xY);
D3Q19_MapRecv(1,-1,0,Dm->recvList_xY,0,recvCount_xY,dvcRecvDist_xY);
//...Map recieve list for the xY edge <<<10)................................
D3Q19_MapRecv(-1,1,0,Dm.recvList_Xy,0,recvCount_Xy,dvcRecvDist_Xy);
D3Q19_MapRecv(-1,1,0,Dm->recvList_Xy,0,recvCount_Xy,dvcRecvDist_Xy);
//...Map recieve list for the XY edge <<<7)................................
D3Q19_MapRecv(1,1,0,Dm.recvList_xy,0,recvCount_xy,dvcRecvDist_xy);
D3Q19_MapRecv(1,1,0,Dm->recvList_xy,0,recvCount_xy,dvcRecvDist_xy);
//...Map recieve list for the xz edge <<<12)................................
D3Q19_MapRecv(-1,0,-1,Dm.recvList_XZ,0,recvCount_XZ,dvcRecvDist_XZ);
D3Q19_MapRecv(-1,0,-1,Dm->recvList_XZ,0,recvCount_XZ,dvcRecvDist_XZ);
//...Map recieve list for the xZ edge <<<14)................................
D3Q19_MapRecv(-1,0,1,Dm.recvList_Xz,0,recvCount_Xz,dvcRecvDist_Xz);
D3Q19_MapRecv(-1,0,1,Dm->recvList_Xz,0,recvCount_Xz,dvcRecvDist_Xz);
//...Map recieve list for the Xz edge <<<13)................................
D3Q19_MapRecv(1,0,-1,Dm.recvList_xZ,0,recvCount_xZ,dvcRecvDist_xZ);
D3Q19_MapRecv(1,0,-1,Dm->recvList_xZ,0,recvCount_xZ,dvcRecvDist_xZ);
//...Map recieve list for the XZ edge <<<11)................................
D3Q19_MapRecv(1,0,1,Dm.recvList_xz,0,recvCount_xz,dvcRecvDist_xz);
D3Q19_MapRecv(1,0,1,Dm->recvList_xz,0,recvCount_xz,dvcRecvDist_xz);
//...Map recieve list for the yz edge <<<16)................................
D3Q19_MapRecv(0,-1,-1,Dm.recvList_YZ,0,recvCount_YZ,dvcRecvDist_YZ);
D3Q19_MapRecv(0,-1,-1,Dm->recvList_YZ,0,recvCount_YZ,dvcRecvDist_YZ);
//...Map recieve list for the yZ edge <<<18)................................
D3Q19_MapRecv(0,-1,1,Dm.recvList_Yz,0,recvCount_Yz,dvcRecvDist_Yz);
D3Q19_MapRecv(0,-1,1,Dm->recvList_Yz,0,recvCount_Yz,dvcRecvDist_Yz);
//...Map recieve list for the Yz edge <<<17)................................
D3Q19_MapRecv(0,1,-1,Dm.recvList_yZ,0,recvCount_yZ,dvcRecvDist_yZ);
D3Q19_MapRecv(0,1,-1,Dm->recvList_yZ,0,recvCount_yZ,dvcRecvDist_yZ);
//...Map recieve list for the YZ edge <<<15)................................
D3Q19_MapRecv(0,1,1,Dm.recvList_yz,0,recvCount_yz,dvcRecvDist_yz);
D3Q19_MapRecv(0,1,1,Dm->recvList_yz,0,recvCount_yz,dvcRecvDist_yz);
//...................................................................................
//......................................................................................

2
common/ScaLBL.h
View File

@@ -129,7 +129,7 @@ extern "C" void ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, int
class ScaLBL_Communicator{
public:
//......................................................................................
ScaLBL_Communicator(Domain &Dm);
ScaLBL_Communicator(Domain *Dm);
//ScaLBL_Communicator(Domain &Dm, IntArray &Map);
~ScaLBL_Communicator();

133
models/ColorModel.cpp
View File

@@ -71,22 +71,17 @@ void ScaLBL_ColorModel::ReadParams(string filename){
MPI_Barrier(comm);
}
ScaLBL_ColorModel::~ScaLBL_ColorModel(){
}
void ScaLBL_ColorModel::ReadInput(){
int rank=Dm->rank();
size_t readID;
//.......................................................................
if (rank == 0) printf("Read input media... \n");
//.......................................................................
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankString,"%05d",Dm->rank());
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
// .......... READ THE INPUT FILE .......................................
id = new char[N];
double sum, sum_local;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
//...........................................................................
if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
//.......................................................................
@@ -98,16 +93,15 @@ void ScaLBL_ColorModel::ReadInput(){
if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
sprintf(LocalRankFilename,"ID.%05i",rank);
size_t readID;
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (IDFILE==NULL) ERROR("lbpm_color_simulator: Error opening file: ID.xxxxx");
readID=fread(id,1,N,IDFILE);
if (readID != size_t(N)) printf("lbpm_color_simulator: Error reading ID (rank=%i) \n",rank);
if (readID != size_t(N)) printf("lbpm_color_simulator: Error reading ID (rank=%i) \n",Dm->rank());
fclose(IDFILE);
// Read restart file
if (Restart == true){
if (rank==0){
if (Dm->rank()==0){
printf("Reading restart file! \n");
ifstream restart("Restart.txt");
if (restart.is_open()){
@@ -136,8 +130,13 @@ void ScaLBL_ColorModel::Create(){
//.......................................................................
// Compute the media porosity, assign phase labels and solid composition
//.......................................................................
double porosity;
double sum_local=0.0;
int nprocs=nprocx*nprocy*nprocz;
int rank=Dm->rank();
double sum;
double porosity;
double sum_local=0.0;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
Np=0; // number of local pore nodes
//.......................................................................
for (int k=1;k<Nz-1;k++){
@@ -204,7 +203,7 @@ void ScaLBL_ColorModel::Create(){
if (rank==0) printf ("Set up memory efficient layout \n");
Map.resize(Nx,Ny,Nz); Map.fill(0);
auto 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);
//...........................................................................
@@ -232,6 +231,9 @@ void ScaLBL_ColorModel::Create(){
//...........................................................................
// Update GPU data structures
if (rank==0) printf ("Setting up device map and neighbor list \n");
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
int *TmpMap;
TmpMap=new int[Np];
for (int k=1; k<Nz-1; k++){
@@ -251,9 +253,12 @@ void ScaLBL_ColorModel::Initialize(){
/*
* This function initializes model (includes both mobile and immobile components)
*/
int rank=Dm->rank();
// Compute the solid interaction potential and copy result to device
if (rank==0) printf("Computing solid interaction potential \n");
double *PhaseLabel;
PhaseLabel=new double [Nx*Ny*Nz];
Mask->AssignComponentLabels(PhaseLabel);
double *Tmp;
Tmp=new double[3*Np];
//Averages->UpdateMeshValues(); // this computes the gradient of distance field (among other things)
@@ -339,9 +344,9 @@ void ScaLBL_ColorModel::Initialize(){
DoubleArray Psy(Nx,Ny,Nz);
DoubleArray Psz(Nx,Ny,Nz);
DoubleArray Psnorm(Nx,Ny,Nz);
ScaLBL_Comm.RegularLayout(Map,&SolidPotential[0],Psx);
ScaLBL_Comm.RegularLayout(Map,&SolidPotential[Np],Psy);
ScaLBL_Comm.RegularLayout(Map,&SolidPotential[2*Np],Psz);
ScaLBL_Comm->RegularLayout(Map,&SolidPotential[0],Psx);
ScaLBL_Comm->RegularLayout(Map,&SolidPotential[Np],Psy);
ScaLBL_Comm->RegularLayout(Map,&SolidPotential[2*Np],Psz);
for (int n=0; n<N; n++) Psnorm(n) = Psx(n)*Psx(n)+Psy(n)*Psy(n)+Psz(n)*Psz(n);
FILE *PFILE;
sprintf(LocalRankFilename,"Potential.%05i.raw",rank);
@@ -367,8 +372,6 @@ void ScaLBL_ColorModel::Initialize(){
}
}
//printf("sw=%f \n",count_wet/double(Np));
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
// initialize phi based on PhaseLabel (include solid component labels)
ScaLBL_CopyToDevice(Phi, PhaseLabel, Np*sizeof(double));
//...........................................................................
@@ -376,11 +379,15 @@ void ScaLBL_ColorModel::Initialize(){
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);
}
void ScaLBL_ColorModel::Run(){
int nprocs=nprocx*nprocy*nprocz;
int rank=Dm->rank();
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
//.......create and start timer............
@@ -390,8 +397,14 @@ void ScaLBL_ColorModel::Run(){
starttime = MPI_Wtime();
//.........................................
//************ MAIN ITERATION LOOP ***************************************/
// read the input database
auto db = std::make_shared<Database>( "input.db" );
auto domain_db = db->getDatabase( "Domain" );
auto color_db = db->getDatabase( "Color" );
auto analysis_db = db->getDatabase( "Analysis" );
PROFILE_START("Loop");
runAnalysis analysis( analysis_db, rank_info, ScaLBL_Comm, Dm, Np, pBC, beta, Map );
// runAnalysis analysis( analysis_db, rank_info, ScaLBL_Comm, Dm, Np, pBC, beta, Map );
while (timestep < timestepMax ) {
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
PROFILE_START("Update");
@@ -399,83 +412,83 @@ void ScaLBL_ColorModel::Run(){
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);
PROFILE_STOP("Update");
// Run the analysis
analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
// analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
}
analysis.finish();
// analysis.finish();
PROFILE_STOP("Loop");
PROFILE_SAVE("lbpm_color_simulator",1);
//************************************************************************

2
tests/CMakeLists.txt
View File

@@ -11,7 +11,7 @@ ADD_LBPM_EXECUTABLE( lbpm_random_pp )
ADD_LBPM_EXECUTABLE( lbpm_refine_pp )
ADD_LBPM_EXECUTABLE( lbpm_morphdrain_pp )
ADD_LBPM_EXECUTABLE( lbpm_morphopen_pp )
ADD_LBPM_EXECUTABLE( lbpm_morph_pp )
#ADD_LBPM_EXECUTABLE( lbpm_morph_pp )
ADD_LBPM_EXECUTABLE( lbpm_segmented_pp )
#ADD_LBPM_EXECUTABLE( lbpm_block_pp )
ADD_LBPM_EXECUTABLE( lbpm_segmented_decomp )