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

eat the cake

Browse Source
This commit is contained in:
James E McClure
2018-05-17 21:03:11 -04:00
parent 84be31aeaf
commit 9e5282fde3
12 changed files with 188 additions and 189 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
tests/CMakeLists.txt
View File

@@ -2,10 +2,10 @@
#ADD_LBPM_EXECUTABLE( lbpm_nonnewtonian_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_nondarcy_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator )
ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )
ADD_LBPM_EXECUTABLE( lbpm_dfh_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )
#ADD_LBPM_EXECUTABLE( lbpm_dfh_simulator )
ADD_LBPM_EXECUTABLE( lbpm_sphere_pp )
ADD_LBPM_EXECUTABLE( lbpm_random_pp )
ADD_LBPM_EXECUTABLE( lbpm_refine_pp )
@@ -41,10 +41,10 @@ ADD_LBPM_TEST( pmmc_cylinder )
ADD_LBPM_TEST( TestTorus )
ADD_LBPM_TEST( TestFluxBC )
ADD_LBPM_TEST( TestMap )
ADD_LBPM_TEST( TestMRT )
ADD_LBPM_TEST( TestColorGrad )
ADD_LBPM_TEST( TestBubbleDFH )
ADD_LBPM_TEST( TestColorGradDFH )
#ADD_LBPM_TEST( TestMRT )
#ADD_LBPM_TEST( TestColorGrad )
#ADD_LBPM_TEST( TestBubbleDFH )
#ADD_LBPM_TEST( TestColorGradDFH )
ADD_LBPM_TEST( TestColorMassBounceback )
ADD_LBPM_TEST( TestPressVel ../example/Piston/input.db)
ADD_LBPM_TEST( TestPoiseuille ../example/Piston/input.db)
@@ -59,7 +59,7 @@ ADD_LBPM_TEST_1_2_4( TestBlobIdentify )
#ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 )
ADD_LBPM_TEST_PARALLEL( TestBlobAnalyze 8 )
ADD_LBPM_TEST_PARALLEL( TestSegDist 8 )
ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 )
#ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 )
#ADD_LBPM_TEST_PARALLEL( TestMassConservationD3Q7 1 )
ADD_LBPM_TEST_1_2_4( testCommunication )
ADD_LBPM_TEST_1_2_4( testUtilities )

36
tests/TestColorSquareTube.cpp
View File

@@ -115,8 +115,8 @@ int main(int argc, char **argv)
double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz;
Domain Dm(db);
std::shared_ptr<Domain> Dm (new Domain(db));
Dm->CommInit(comm);
Nx += 2;
Ny += 2;
@@ -136,7 +136,7 @@ int main(int argc, char **argv)
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nx + i;
Dm.id[n]=0;
Dm->id[n]=0;
}
}
}
@@ -155,16 +155,16 @@ int main(int argc, char **argv)
int kglobal= k+(Nz-2)*kproc;
// Initialize phase position field for parallel bubble test
if (iglobal < 2) Dm.id[n]=0;
else if (iglobal > (Nx-2)*nprocx-2) Dm.id[n]=0;
else if (jglobal < 2) Dm.id[n]=0;
else if (jglobal > (Ny-2)*nprocy-2) Dm.id[n]=0;
else if (kglobal < 20) Dm.id[n]=1;
else Dm.id[n]=2;
if (iglobal < 2) Dm->id[n]=0;
else if (iglobal > (Nx-2)*nprocx-2) Dm->id[n]=0;
else if (jglobal < 2) Dm->id[n]=0;
else if (jglobal > (Ny-2)*nprocy-2) Dm->id[n]=0;
else if (kglobal < 20) Dm->id[n]=1;
else Dm->id[n]=2;
}
}
}
Dm.CommInit(comm);
Dm->CommInit(comm);
//.......................................................................
// Compute the media porosity, assign phase labels and solid composition
@@ -174,13 +174,13 @@ int main(int argc, char **argv)
int Np=0; // number of local pore nodes
double *PhaseLabel;
PhaseLabel = new double[N];
Dm.AssignComponentLabels(PhaseLabel);
Dm->AssignComponentLabels(PhaseLabel);
//.......................................................................
for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){
for (i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (Dm.id[n] > 0){
if (Dm->id[n] > 0){
sum_local+=1.0;
Np++;
}
@@ -203,7 +203,7 @@ int main(int argc, char **argv)
char *ID;
ScaLBL_AllocateDeviceMemory((void **) &ID, N); // Allocate device memory
// Copy to the device
ScaLBL_CopyToDevice(ID, Dm.id, N);
ScaLBL_CopyToDevice(ID, Dm->id, N);
//...........................................................................
if (rank==0){
@@ -219,7 +219,7 @@ int main(int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Np];
ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np);
MPI_Barrier(comm);
//......................device distributions.................................
@@ -277,12 +277,12 @@ int main(int argc, char **argv)
if (rank==0) printf ("Initializing distributions \n");
// Initialize the phase field and variables
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm.last_interior, Np);
if (Dm.kproc()==0){
if (Dm->kproc()==0){
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
}
if (Dm.kproc() == nprocz-1){
if (Dm->kproc() == nprocz-1){
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
@@ -409,7 +409,7 @@ 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;
if (Dm.id[n] > 0){
if (Dm->id[n] > 0){
int idx = Map(i,j,k);
sum_local+=VEL[2*Np+idx];
}
@@ -426,7 +426,7 @@ 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;
if (Dm.id[n] > 0){
if (Dm->id[n] > 0){
int idx = Map(i,j,k);
double vz = VEL[2*Np+idx];
printf("%f ",vz);

46
tests/TestFluxBC.cpp
View File

@@ -31,7 +31,7 @@ int main (int argc, char **argv)
double din,dout;
int BC=1;
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));
Nz += 2;
Nx = Ny = Nz; // Cubic domain
@@ -67,8 +67,8 @@ int main (int argc, char **argv)
id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
//.........................................................
// Initialize communication structures in averaging domain
for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = id[i];
Dm.CommInit(comm);
for (i=0; i<Dm->Nx*Dm->Ny*Dm->Nz; i++) Dm->id[i] = id[i];
Dm->CommInit(comm);
Np=0; // number of local pore nodes
//.......................................................................
@@ -76,7 +76,7 @@ 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;
if (Dm.id[n] > 0){
if (Dm->id[n] > 0){
Np++;
}
}
@@ -85,7 +85,7 @@ int main (int argc, char **argv)
//...........................................................................
if (rank==0) printf ("Create ScaLBL_Communicator \n");
// Create a communicator for the device
ScaLBL_Communicator ScaLBL_Comm(Dm);
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm(new ScaLBL_Communicator(Dm));
if (rank==0) printf ("Set up the neighborlist \n");
int Npad=Np+32;
@@ -94,7 +94,7 @@ int main (int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np);
MPI_Barrier(comm);
//......................device distributions.................................
@@ -142,7 +142,7 @@ int main (int argc, char **argv)
double flux = 1.0;
int timestep=0;
din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
printf("Computed pressure for flux = %f\n",din);
@@ -165,7 +165,7 @@ 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;
if (Dm.id[n] > 0){
if (Dm->id[n] > 0){
int idx = Map(i,j,k);
Q += VEL[2*Np+idx];
}
@@ -194,21 +194,21 @@ int main (int argc, char **argv)
while (timestep < 2000) {
ScaLBL_Comm.SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, 0, ScaLBL_Comm->next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
ScaLBL_Comm.SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_MRT(fq, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
ScaLBL_D3Q19_AAeven_MRT(fq, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_MRT(fq, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
ScaLBL_D3Q19_AAeven_MRT(fq, 0, ScaLBL_Comm->next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
//************************************************************************/
@@ -226,7 +226,7 @@ int main (int argc, char **argv)
for (k=1;k<Nz-1;k++){
for (i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
if (Dm.id[n] > 0){
if (Dm->id[n] > 0){
int idx = Map(i,j,k);
double vz = VEL[2*Np+idx];
printf("%f ",vz);
@@ -242,7 +242,7 @@ 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;
if (Dm.id[n] > 0){
if (Dm->id[n] > 0){
int idx = Map(i,j,k);
double vz = VEL[2*Np+idx];
printf("%f ",vz);
@@ -256,7 +256,7 @@ 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;
if (Dm.id[n] > 0){
if (Dm->id[n] > 0){
int idx = Map(i,j,k);
Q += VEL[2*Np+idx];
}

49
tests/TestForceMoments.cpp
View File

@@ -13,7 +13,6 @@ using namespace std;
extern void PrintNeighborList(int * neighborList, int Np, int rank) {
if (rank == 0) {
int n;
int neighbor;
for (int i = 0; i < Np; i++) {
@@ -32,7 +31,7 @@ extern void PrintNeighborList(int * neighborList, int Np, int rank) {
std::shared_ptr<Database> loadInputs( int nprocs )
{
auto db = std::make_shared<Database>( "Domain.in" );
const int dim = 50;
db->putScalar<int>( "BC", 0 );
if ( nprocs == 1 ){
db->putVector<int>( "nproc", { 1, 1, 1 } );
@@ -118,7 +117,7 @@ int main(int argc, char **argv)
double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz;
Domain Dm(domain_db);
std::shared_ptr<Domain> Dm(new Domain(domain_db));
Nx += 2;
Ny += 2;
@@ -135,7 +134,7 @@ int main(int argc, char **argv)
/*
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx");
fread(Dm.id,1,N,IDFILE);
fread(Dm->id,1,N,IDFILE);
fclose(IDFILE);
*/
@@ -144,11 +143,11 @@ int main(int argc, char **argv)
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny+j*Nx+i;
Dm.id[n]=1;
Dm->id[n]=1;
}
}
}
Dm.CommInit(comm);
Dm->CommInit(comm);
MPI_Barrier(comm);
if (rank == 0) cout << "Domain set." << endl;
@@ -157,7 +156,7 @@ 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;
if (Dm.id[n] > 0){
if (Dm->id[n] > 0){
Np++;
}
}
@@ -165,16 +164,14 @@ int main(int argc, char **argv)
}
if (rank==0) printf ("Create ScaLBL_Communicator \n");
// Create a communicator for the device
ScaLBL_Communicator ScaLBL_Comm(Dm);
auto ScaLBL_Comm = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Dm));
//...........device phase ID.................................................
if (rank==0) printf ("Copying phase ID to device \n");
char *ID;
ScaLBL_AllocateDeviceMemory((void **) &ID, N); // Allocate device memory
// Copy to the device
ScaLBL_CopyToDevice(ID, Dm.id, N);
ScaLBL_CopyToDevice(ID, Dm->id, N);
//...........................................................................
if (rank==0){
@@ -191,7 +188,7 @@ int main(int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Np];
ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np);
if (rank == 0) PrintNeighborList(neighborList,Np, rank);
@@ -225,8 +222,8 @@ int main(int argc, char **argv)
starttime = MPI_Wtime();
/************ MAIN ITERATION LOOP (timing communications)***************************************/
//ScaLBL_Comm.SendD3Q19(dist, &dist[10*Np]);
//ScaLBL_Comm.RecvD3Q19(dist, &dist[10*Np]);
//ScaLBL_Comm->SendD3Q19(dist, &dist[10*Np]);
//ScaLBL_Comm->RecvD3Q19(dist, &dist[10*Np]);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
if (rank==0) printf("Beginning AA timesteps...\n");
@@ -235,17 +232,17 @@ int main(int argc, char **argv)
while (timestep < 2) {
ScaLBL_Comm.SendD3Q19AA(dist); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm.RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->SendD3Q19AA(dist); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, 0, ScaLBL_Comm->next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
ScaLBL_Comm.SendD3Q19AA(dist); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_MRT(dist, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm.RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAeven_MRT(dist, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->SendD3Q19AA(dist); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_MRT(dist, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_Comm->RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
ScaLBL_D3Q19_AAeven_MRT(dist, 0, ScaLBL_Comm->next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
timestep++;
//************************************************************************/
@@ -293,7 +290,7 @@ 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;
//printf("%i ",Dm.id[n]);
//printf("%i ",Dm->id[n]);
n = Map(i,j,k);
//printf("%i,%i,%i; %i :",i,j,k,n);
if (n<0) vz =0.f;
@@ -318,7 +315,7 @@ int main(int argc, char **argv)
for (k=1;k<Nz-1;k++){
for (j=1;j<Ny-1;j++){
n = k*Nx*Ny+j*Nx+i;
//printf("%i ",Dm.id[n]);
//printf("%i ",Dm->id[n]);
n = Map(i,j,k);
double fa = DIST[(2*q+1)*Np+n];
double fb = DIST[2*(q+1)*Np+n];
@@ -339,7 +336,7 @@ int main(int argc, char **argv)
for (k=1;k<Nz-1;k++){
for (i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
//printf("%i ",Dm.id[n]);
//printf("%i ",Dm->id[n]);
n = Map(i,j,k);
double fa = DIST[(2*q+1)*Np+n];
double fb = DIST[2*(q+1)*Np+n];
@@ -359,7 +356,7 @@ 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;
//printf("%i ",Dm.id[n]);
//printf("%i ",Dm->id[n]);
n = Map(i,j,k);
double fa = DIST[(2*q+1)*Np+n];
double fb = DIST[2*(q+1)*Np+n];

15
tests/TestMap.cpp
View File

@@ -134,7 +134,7 @@ int main(int argc, char **argv)
MPI_Barrier(comm);
int BoundaryCondition=0;
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
std::shared_ptr<Domain> Dm(new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition));
Nx += 2;
Ny += 2;
@@ -147,18 +147,17 @@ 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;
Np++;
}
}
}
Dm.CommInit(comm);
Dm->CommInit(comm);
// Create a communicator for the device (will use optimized layout)
ScaLBL_Communicator ScaLBL_Comm(Dm);
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm(new ScaLBL_Communicator(Dm));
//Create a second communicator based on the regular data layout
ScaLBL_Communicator ScaLBL_Comm_Regular(Dm);
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm_Regular(new ScaLBL_Communicator(Dm));
if (rank==0){
printf("Total domain size = %i \n",N);
@@ -173,7 +172,7 @@ int main(int argc, char **argv)
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np);
MPI_Barrier(comm);
//......................device distributions.................................
@@ -236,7 +235,7 @@ int main(int argc, char **argv)
// Loop over the distributions for interior lattice sites
if (rank==0) printf ("Loop over distributions \n");
for (int idx=ScaLBL_Comm.first_interior; idx<ScaLBL_Comm.last_interior; idx++){
for (int idx=ScaLBL_Comm->first_interior; idx<ScaLBL_Comm->last_interior; idx++){
n = TmpMap[idx];
k = n/(Nx*Ny);
j = (n-Nx*Ny*k)/Nx;

16
tests/TestPressVel.cpp
View File

@@ -76,7 +76,8 @@ int main(int argc, char **argv)
double iVol_global = 1.0/Nx/Ny/Nz/nprocx/nprocy/nprocz;
Domain Dm(domain_db);
std::shared_ptr<Domain> Dm (new Domain(domain_db));
Dm->CommInit(comm);
Nx += 2;
Ny += 2;
@@ -91,7 +92,7 @@ int main(int argc, char **argv)
char LocalRankFilename[40];
sprintf(LocalRankFilename,"ID.%05i",rank);
Dm.CommInit(comm);
Dm->CommInit(comm);
//.......................................................................
// Compute the media porosity
@@ -103,8 +104,8 @@ 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;
if (Dm.id[n] > 0){
Dm->id[n] = 1;
if (Dm->id[n] > 0){
sum_local+=1.0;
Np++;
}
@@ -120,14 +121,13 @@ int main(int argc, char **argv)
if (rank==0) printf ("Create ScaLBL_Communicator \n");
// Create a communicator for the device
ScaLBL_Communicator ScaLBL_Comm(Dm);
auto ScaLBL_Comm = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Dm));
//...........device phase ID.................................................
if (rank==0) printf ("Copying phase ID to device \n");
char *ID;
ScaLBL_AllocateDeviceMemory((void **) &ID, N); // Allocate device memory
// Copy to the device
ScaLBL_CopyToDevice(ID, Dm.id, N);
ScaLBL_CopyToDevice(ID, Dm->id, N);
//...........................................................................
if (rank==0){
@@ -142,7 +142,7 @@ int main(int argc, char **argv)
int *neighborList;
IntArray Map(Nx,Ny,Nz);
neighborList= new int[18*Npad];
Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Dm->id,Np);
MPI_Barrier(comm);
//......................device distributions.................................

82
tests/TestTorus.cpp
View File

@@ -61,9 +61,9 @@ int main(int argc, char **argv)
// Get the rank info
Domain Dm(db);
std::shared_ptr<Domain> Dm(new Domain(db));
// const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
TwoPhase Averages(Dm);
std::shared_ptr<TwoPhase> Averages(new TwoPhase(Dm));
Nx += 2;
Ny += 2;
Nz += 2;
@@ -73,12 +73,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
//.......................................................................
//.......................................................................
@@ -103,42 +103,42 @@ int main(int argc, char **argv)
n = k*Nx*Ny+j*Nx+i;
// global position relative to center
x = Dm.iproc()*Nx+i - CX;
y = Dm.jproc()*Ny+j - CY;
z = Dm.kproc()*Nz+k - CZ;
x = Dm->iproc()*Nx+i - CX;
y = Dm->jproc()*Ny+j - CY;
z = Dm->kproc()*Nz+k - CZ;
// Shrink the sphere sizes by two voxels to make sure they don't touch
Averages.SDs(i,j,k) = 100.0;
Averages->SDs(i,j,k) = 100.0;
//..............................................................................
// Single torus
Averages.Phase(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2;
Averages->Phase(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2;
// Double torus
/* y = Dm.jproc()*Ny+j - CY1;
//z = Dm.kproc()*Nz+k - CZ +R1;
Averages.Phase(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2;
/* y = Dm->jproc()*Ny+j - CY1;
//z = Dm->kproc()*Nz+k - CZ +R1;
Averages->Phase(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2;
y = Dm.jproc()*Ny+j - CY2;
//z = Dm.kproc()*Nz+k - CZ-R1;
Averages.Phase(i,j,k) = min(Averages.Phase(i,j,k),
y = Dm->jproc()*Ny+j - CY2;
//z = Dm->kproc()*Nz+k - CZ-R1;
Averages->Phase(i,j,k) = min(Averages->Phase(i,j,k),
sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2);
*///..............................................................................
//Averages.Phase(i,j,k) = - Averages.Phase(i,j,k);
if (Averages.Phase(i,j,k) > 0.0){
Dm.id[n] = 2;
//Averages->Phase(i,j,k) = - Averages->Phase(i,j,k);
if (Averages->Phase(i,j,k) > 0.0){
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);
Averages.Phase_tplus(i,j,k) = Averages.SDn(i,j,k);
Averages.Phase_tminus(i,j,k) = Averages.SDn(i,j,k);
Averages.DelPhi(i,j,k) = 0.0;
Averages.Press(i,j,k) = 0.0;
Averages.Vel_x(i,j,k) = 0.0;
Averages.Vel_y(i,j,k) = 0.0;
Averages.Vel_z(i,j,k) = 0.0;
Averages->SDn(i,j,k) = Averages->Phase(i,j,k);
Averages->Phase(i,j,k) = Averages->SDn(i,j,k);
Averages->Phase_tplus(i,j,k) = Averages->SDn(i,j,k);
Averages->Phase_tminus(i,j,k) = Averages->SDn(i,j,k);
Averages->DelPhi(i,j,k) = 0.0;
Averages->Press(i,j,k) = 0.0;
Averages->Vel_x(i,j,k) = 0.0;
Averages->Vel_y(i,j,k) = 0.0;
Averages->Vel_z(i,j,k) = 0.0;
}
}
}
@@ -146,25 +146,25 @@ int main(int argc, char **argv)
double beta = 0.95;
if (rank==0) printf("initializing the system \n");
Averages.UpdateSolid();
Dm.CommunicateMeshHalo(Averages.Phase);
Dm.CommunicateMeshHalo(Averages.SDn);
Averages->UpdateSolid();
Dm->CommunicateMeshHalo(Averages->Phase);
Dm->CommunicateMeshHalo(Averages->SDn);
Averages.Initialize();
Averages.UpdateMeshValues();
Averages->Initialize();
Averages->UpdateMeshValues();
if (rank==0) printf("computing local averages \n");
Averages.AssignComponentLabels();
Averages.ComponentAverages();
Averages.PrintComponents(int(5));
Averages->AssignComponentLabels();
Averages->ComponentAverages();
Averages->PrintComponents(int(5));
if (rank==0) printf("reducing averages \n");
Averages.Initialize();
Averages.ComputeLocal();
Averages.Reduce();
Averages.PrintAll(int(5));
// Averages.Reduce();
Averages->Initialize();
Averages->ComputeLocal();
Averages->Reduce();
Averages->PrintAll(int(5));
// Averages->Reduce();
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Barrier(comm);

24
tests/lbpm_captube_pp.cpp
View File

@@ -34,8 +34,6 @@ int main(int argc, char **argv)
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
{
// parallel domain size (# of sub-domains)
int sendtag,recvtag;
//*****************************************
// MPI ranks for all 18 neighbors
//**********************************
@@ -44,8 +42,6 @@ int main(int argc, char **argv)
int rank_xz,rank_XZ,rank_xZ,rank_Xz;
int rank_yz,rank_YZ,rank_yZ,rank_Yz;
//**********************************
MPI_Request req1[18],req2[18];
MPI_Status stat1[18],stat2[18];
double TubeRadius =15.0;
int BC;
@@ -88,10 +84,12 @@ int main(int argc, char **argv)
printf("********************************************************\n");
}
// Initialized domain and averaging framework for Two-Phase Flow
Domain Dm(db);
Dm.CommInit(comm);
TwoPhase Averages(Dm);
double Lx=1.f;
double Ly=1.f;
double Lz=1.f;
std::shared_ptr<Domain> Dm (new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
Dm->CommInit(comm);
std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc,
rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z,
@@ -133,18 +131,18 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
// Cylindrical capillary tube aligned with the z direction
Averages.SDs(i,j,k) = TubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)
Averages->SDs(i,j,k) = TubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)
+ (j-Ny/2)*(j-Ny/2)));
// Initialize phase positions
if (Averages.SDs(i,j,k) < 0.0){
if (Averages->SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.kproc()*Nz+k<BubbleBottom){
else if (Dm->kproc()*Nz+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.kproc()*Nz+k<BubbleTop){
else if (Dm->kproc()*Nz+k<BubbleTop){
id[n] = 1;
sum++;
}
@@ -177,7 +175,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fwrite(Averages->SDs.data(),8,Averages->SDs.length(),DIST);
fclose(DIST);
sprintf(LocalRankFilename,"ID.%05i",rank);

23
tests/lbpm_inkbottle_pp.cpp
View File

@@ -26,7 +26,6 @@ int main(int argc, char **argv)
// parallel domain size (# of sub-domains)
int nprocx,nprocy,nprocz;
int iproc,jproc,kproc;
int sendtag,recvtag;
//*****************************************
// MPI ranks for all 18 neighbors
//**********************************
@@ -35,8 +34,6 @@ int main(int argc, char **argv)
int rank_xz,rank_XZ,rank_xZ,rank_Xz;
int rank_yz,rank_YZ,rank_yZ,rank_Yz;
//**********************************
MPI_Request req1[18],req2[18];
MPI_Status stat1[18],stat2[18];
double UpperTubeRadius =15.0;
double LowerTubeRadius =15.0;
@@ -117,9 +114,9 @@ int main(int argc, char **argv)
}
// Initialized domain and averaging framework for Two-Phase Flow
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
Dm.CommInit(comm);
TwoPhase Averages(Dm);
std::shared_ptr<Domain> Dm (new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
Dm->CommInit(comm);
std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc,
rank_x, rank_y, rank_z, rank_X, rank_Y, rank_Z,
@@ -163,25 +160,25 @@ int main(int argc, char **argv)
n = k*Nx*Ny + j*Nz + i;
// Cylindrical capillary tube aligned with the z direction
if (k<Nz/2)
Averages.SDs(i,j,k) = LowerTubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)
Averages->SDs(i,j,k) = LowerTubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)
+ (j-Ny/2)*(j-Ny/2)));
else
Averages.SDs(i,j,k) = UpperTubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)
Averages->SDs(i,j,k) = UpperTubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)
+ (j-Ny/2)*(j-Ny/2)));
BulbDist = BulbRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)+ (j-Ny/2)*(j-Ny/2) + (k-Nz/2)*(k-Nz/2)));
if (BulbDist > Averages.SDs(i,j,k)) Averages.SDs(i,j,k) = BulbDist;
if (BulbDist > Averages->SDs(i,j,k)) Averages->SDs(i,j,k) = BulbDist;
// Initialize phase positions
if (Averages.SDs(i,j,k) < 0.0){
if (Averages->SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Dm.kproc()*Nz+k<BubbleBottom){
else if (Dm->kproc()*Nz+k<BubbleBottom){
id[n] = 2;
sum++;
}
else if (Dm.kproc()*Nz+k<BubbleTop){
else if (Dm->kproc()*Nz+k<BubbleTop){
id[n] = 1;
sum++;
}
@@ -214,7 +211,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fwrite(Averages->SDs.data(),8,Averages->SDs.length(),DIST);
fclose(DIST);
sprintf(LocalRankFilename,"ID.%05i",rank);

16
tests/lbpm_plates_pp.cpp
View File

@@ -110,9 +110,11 @@ int main(int argc, char **argv)
}
// Initialized domain and averaging framework for Two-Phase Flow
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
Dm.CommInit(comm);
TwoPhase Averages(Dm);
std::shared_ptr<Domain> Dm (new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
Dm->CommInit(comm);
std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
MPI_Barrier(comm);
@@ -149,13 +151,13 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
// Cylindrical capillary tube aligned with the z direction
Averages.SDs(i,j,k) = TubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)));
Averages->SDs(i,j,k) = TubeRadius-sqrt(1.0*((i-Nx/2)*(i-Nx/2)));
// Initialize phase positions
if (Averages.SDs(i,j,k) < 0.0){
if (Averages->SDs(i,j,k) < 0.0){
id[n] = 0;
}
else if (Averages.SDs(i,j,k) < WIDTH){
else if (Averages->SDs(i,j,k) < WIDTH){
id[n] = 2;
sum++;
}
@@ -188,7 +190,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fwrite(Averages->SDs.data(),8,Averages->SDs.length(),DIST);
fclose(DIST);
sprintf(LocalRankFilename,"ID.%05i",rank);

22
tests/lbpm_porenetwork_pp.cpp
View File

@@ -100,9 +100,13 @@ int main(int argc, char **argv)
}
// Initialized domain and averaging framework for Two-Phase Flow
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
Dm.CommInit(comm);
TwoPhase Averages(Dm);
// Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
//Dm.CommInit(comm);
//TwoPhase Averages(Dm);
std::shared_ptr<Domain> Dm (new Domain(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
Dm->CommInit(comm);
std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
MPI_Barrier(comm);
@@ -180,7 +184,7 @@ int main(int argc, char **argv)
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
Averages.SDs(i,j,k) = -2.f*Nx;
Averages->SDs(i,j,k) = -2.f*Nx;
}
}
}
@@ -208,7 +212,7 @@ int main(int argc, char **argv)
double zk = double(k) - z;
// value of s along center line {x=alpha*s, y = beta*s, z=gamma*s} that is closest to i,j,k
double s = (alpha*xi + beta*yj + gamma*zk)/(alpha*alpha + beta*beta + gamma*gamma);
double distance=Averages.SDs(i,j,k);
double distance=Averages->SDs(i,j,k);
if (s > length){
//distance = radius - sqrt((xi-X)*(xi-X) + (yj-Y)*(yj-Y) + (zk-Z)*(zk-Z));
}
@@ -222,7 +226,7 @@ int main(int argc, char **argv)
distance = radius - sqrt((xi-xs)*(xi-xs) + (yj-ys)*(yj-ys) + (zk-zs)*(zk-zs));
//if (distance>0)printf("s=%f,alpha=%f,beta=%f,gamma=%f,distance=%f\n",s,alpha,beta,gamma,distance);
}
if (distance > Averages.SDs(i,j,k)) Averages.SDs(i,j,k) = distance;
if (distance > Averages->SDs(i,j,k)) Averages->SDs(i,j,k) = distance;
}
// Compute the distance to each sphere
@@ -237,7 +241,7 @@ int main(int argc, char **argv)
double distance = radius - sqrt((xi-x)*(xi-x) + (yj-y)*(yj-y) + (zk-z)*(zk-z));
if (distance > Averages.SDs(i,j,k)) Averages.SDs(i,j,k) = distance;
if (distance > Averages->SDs(i,j,k)) Averages->SDs(i,j,k) = distance;
}
}
}
@@ -247,7 +251,7 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
// Initialize phase positions
if (Averages.SDs(i,j,k) < 0.0){
if (Averages->SDs(i,j,k) < 0.0){
id[n] = 0;
}
else{
@@ -279,7 +283,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fwrite(Averages->SDs.data(),8,Averages->SDs.length(),DIST);
fclose(DIST);
sprintf(LocalRankFilename,"ID.%05i",rank);

30
tests/lbpm_segmented_pp.cpp
View File

@@ -189,9 +189,11 @@ int main(int argc, char **argv)
char LocalRankFilename[40];
int N = (nx+2)*(ny+2)*(nz+2);
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
for (n=0; n<N; n++) Dm.id[n]=1;
Dm.CommInit(comm);
std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
for (n=0; n<N; n++) Dm->id[n]=1;
Dm->CommInit(comm);
std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
nx+=2; ny+=2; nz+=2;
// Read the phase ID
@@ -200,7 +202,7 @@ int main(int argc, char **argv)
if (READ_FROM_BLOCK == 0){
size_t readID;
FILE *ID = fopen(LocalRankFilename,"rb");
readID=fread(Dm.id,1,N,ID);
readID=fread(Dm->id,1,N,ID);
if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID \n");
fclose(ID);
}
@@ -208,13 +210,13 @@ int main(int argc, char **argv)
// Read from the large block and write the local ID file
if (rank==0) printf("Reading ID file from blocks \n");
fflush(stdout);
porosity = ReadFromBlock(Dm.id,Dm.iproc(),Dm.jproc(),Dm.kproc(),nx,ny,nz);
porosity = ReadFromBlock(Dm->id,Dm->iproc(),Dm->jproc(),Dm->kproc(),nx,ny,nz);
MPI_Barrier(MPI_COMM_WORLD);
if (rank==0) printf("Writing local ID files (poros=%f) \n",porosity);
fflush(stdout);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(Dm.id,1,N,ID);
fwrite(Dm->id,1,N,ID);
fclose(ID);
if (rank==0) printf("Succeeded! \n");
fflush(stdout);
@@ -222,8 +224,8 @@ int main(int argc, char **argv)
}
// Initialize the domain and communication
Array<char> id(nx,ny,nz);
TwoPhase Averages(Dm);
Array<char> id(nx,ny,nz);
//TwoPhase Averages(Dm);
// DoubleArray Distance(nx,ny,nz);
// DoubleArray Phase(nx,ny,nz);
@@ -234,7 +236,7 @@ int main(int argc, char **argv)
for (i=0;i<nx;i++){
n = k*nx*ny+j*nx+i;
// Initialize the solid phase
if (Dm.id[n] == 0) id(i,j,k) = 0;
if (Dm->id[n] == 0) id(i,j,k) = 0;
else id(i,j,k) = 1;
}
}
@@ -245,17 +247,17 @@ int main(int argc, char **argv)
for (i=0;i<nx;i++){
n=k*nx*ny+j*nx+i;
// Initialize distance to +/- 1
Averages.SDs(i,j,k) = 2.0*double(id(i,j,k))-1.0;
Averages->SDs(i,j,k) = 2.0*double(id(i,j,k))-1.0;
}
}
}
MeanFilter(Averages.SDs);
MeanFilter(Averages->SDs);
double LocalVar, TotalVar;
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
int Maxtime=10*max(max(Dm.Nx*Dm.nprocx(),Dm.Ny*Dm.nprocy()),Dm.Nz*Dm.nprocz());
int Maxtime=10*max(max(Dm->Nx*Dm->nprocx(),Dm->Ny*Dm->nprocy()),Dm->Nz*Dm->nprocz());
Maxtime=min(Maxtime,MAXTIME);
LocalVar = Eikonal(Averages.SDs,id,Dm,Maxtime);
LocalVar = Eikonal(Averages->SDs,id,*Dm,Maxtime);
MPI_Allreduce(&LocalVar,&TotalVar,1,MPI_DOUBLE,MPI_SUM,comm);
TotalVar /= nprocs;
@@ -263,7 +265,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.data(),8,N,DIST);
fwrite(Averages->SDs.data(),8,N,DIST);
fclose(DIST);
}