Fixed merge conflict in common/TwoPhase.cpp

This commit is contained in:
James E McClure 2015-09-02 08:22:32 -04:00
commit 8178f2b776
5 changed files with 73 additions and 38 deletions

View File

@ -9,7 +9,7 @@
#define GET_ARRAY_INDEX(i1,i2,i3,i4) i1+d_N[0]*(i2+d_N[1]*(i3+d_N[2]*i4))
#ifdef DEBUG
#if defined(DEBUG) || defined(_DEBUG)
#define CHECK_ARRAY_INDEX(i1,i2,i3,i4) \
if ( GET_ARRAY_INDEX(i1,i2,i3,i4)>d_length ) \
ERROR("Index exceeds array bounds");

View File

@ -59,10 +59,10 @@ TwoPhase::TwoPhase(Domain &dm):
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),
TIMELOG(NULL), NWPLOG(NULL), WPLOG(NULL),
Dm(dm), NumberComponents_WP(0), NumberComponents_NWP(0), trimdist(0),
Dm(dm), NumberComponents_WP(0), NumberComponents_NWP(0), trimdist(0),
porosity(0), poreVol(0), awn(0), ans(0), aws(0), lwns(0), wp_volume(0), nwp_volume(0),
As(0), dummy(0), vol_w(0), vol_n(0), sat_w(0), sat_w_previous(0),
pan(0) ,paw(0), pan_global(0), paw_global(0), vol_w_global(0), vol_n_global(0),
pan(0), paw(0), pan_global(0), paw_global(0), vol_w_global(0), vol_n_global(0),
awn_global(0), ans_global(0),aws_global(0), lwns_global(0), efawns(0), efawns_global(0),
Jwn(0), Jwn_global(0), Kwn(0), Kwn_global(0), KNwns(0), KNwns_global(0),
KGwns(0), KGwns_global(0), trawn(0), trawn_global(0), trJwn(0), trJwn_global(0),
@ -366,6 +366,14 @@ void TwoPhase::ComputeLocal()
if (Dm.BoundaryCondition > 0 && Dm.kproc == 0) kmin=4;
if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4;
// Map solid to erode the fluid so that interfaces can be calculated accurately
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
SDs(i,j,k) += 1.0;
}
}
}
for (k=kmin; k<kmax; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
@ -464,6 +472,15 @@ void TwoPhase::ComputeLocal()
}
}
}
// Map solid back
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
SDs(i,j,k) -= 1.0;
}
}
}
}
@ -516,6 +533,15 @@ void TwoPhase::ComponentAverages()
if (Dm.BoundaryCondition > 0 && Dm.kproc == 0) kmin=4;
if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4;
// Map solid to erode the fluid so that interfaces can be calculated accurately
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
SDs(i,j,k) += 1.0;
}
}
}
for (k=kmin; k<kmax; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
@ -712,6 +738,15 @@ void TwoPhase::ComponentAverages()
}
}
}
// Map solid to erode the fluid so that interfaces can be calculated accurately
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
SDs(i,j,k) -= 1.0;
}
}
}
if (Dm.rank==0){
printf("Component averages computed locally -- reducing result... \n");
@ -728,8 +763,7 @@ void TwoPhase::ComponentAverages()
MPI_Barrier(Dm.Comm);
// MPI_Allreduce(&ComponentAverages_NWP(0,0),&RecvBuffer(0,0),BLOB_AVG_COUNT*NumberComponents_NWP,
// MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Reduce(&ComponentAverages_NWP(0,0),&RecvBuffer(0),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
MPI_Reduce(ComponentAverages_NWP.get(),RecvBuffer.get(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
if (Dm.rank==0){
printf("rescaling... \n");
@ -824,9 +858,8 @@ void TwoPhase::ComponentAverages()
// reduce the wetting phase averages
for (int b=0; b<NumberComponents_WP; b++){
MPI_Barrier(Dm.Comm);
// MPI_Allreduce(&ComponentAverages_WP(0,b),&RecvBuffer(0),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Reduce(&ComponentAverages_WP(0,b),&RecvBuffer(0),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
// MPI_Allreduce(&ComponentAverages_WP(0,b),RecvBuffer.get(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Reduce(&ComponentAverages_WP(0,b),RecvBuffer.get(),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);
}

View File

@ -6,6 +6,7 @@
extern "C" void AllocateDeviceMemory(void** address, size_t size){
//cudaMalloc(address,size);
(*address) = malloc(size);
memset(*address,0,size);
if (*address==NULL){
printf("Memory allocation failed! \n");

View File

@ -3,6 +3,7 @@
extern "C" void dvc_AllocateDeviceMemory(void** address, size_t size){
cudaMalloc(address,size);
cudaMemset(*address,0,size);
}
extern "C" void dvc_CopyToDevice(void* dest, void* source, size_t size){

View File

@ -303,7 +303,7 @@ int main(int argc, char **argv)
NULL_USE(pBC); NULL_USE(velBC);
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
TwoPhase Averages(Dm);
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,
@ -350,7 +350,7 @@ int main(int argc, char **argv)
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// WriteLocalSolidID(LocalRankFilename, id, N);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Averages.SDs.get(), N);
ReadBinaryFile(LocalRankFilename, Averages->SDs.get(), N);
MPI_Barrier(MPI_COMM_WORLD);
if (rank == 0) cout << "Domain set." << endl;
@ -371,11 +371,11 @@ int main(int argc, char **argv)
for ( j=0;j<Ny;j++){
for ( i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
if (Averages.SDs(n) > 0.0){
if (Averages->SDs(n) > 0.0){
id[n] = 2;
}
// compute the porosity (actual interface location used)
if (Averages.SDs(n) > 0.0){
if (Averages->SDs(n) > 0.0){
sum++;
}
}
@ -435,7 +435,7 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
//id[n] = 1;
Averages.SDs(n) = max(Averages.SDs(n),1.0*(2.5-k));
Averages->SDs(n) = max(Averages->SDs(n),1.0*(2.5-k));
}
}
}
@ -446,7 +446,7 @@ int main(int argc, char **argv)
for (i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
//id[n] = 2;
Averages.SDs(n) = max(Averages.SDs(n),1.0*(k-Nz+2.5));
Averages->SDs(n) = max(Averages->SDs(n),1.0*(k-Nz+2.5));
}
}
}
@ -528,7 +528,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy signed distance for device initialization
CopyToDevice(dvcSignDist, Averages.SDs.get(), dist_mem_size);
CopyToDevice(dvcSignDist, Averages->SDs.get(), dist_mem_size);
//...........................................................................
int logcount = 0; // number of surface write-outs
@ -571,11 +571,11 @@ int main(int argc, char **argv)
MPI_Barrier(MPI_COMM_WORLD);
//.......................................................................
// Once phase has been initialized, map solid to account for 'smeared' interface
for (i=0; i<N; i++) Averages.SDs(i) -= (1.0); //
//for (i=0; i<N; i++) Averages->SDs(i) -= (1.0); //
//.......................................................................
// Finalize setup for averaging domain
//Averages.SetupCubes(Dm);
Averages.UpdateSolid();
//Averages->SetupCubes(Dm);
Averages->UpdateSolid();
//.......................................................................
//*************************************************************************
@ -638,7 +638,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Averages.Phase_tplus.get(),Phi,N*sizeof(double));
CopyToHost(Averages->Phase_tplus.get(),Phi,N*sizeof(double));
//...........................................................................
//...........................................................................
// Copy the data for for the analysis timestep
@ -647,11 +647,11 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
CopyToHost(Averages->Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages->Press.get(),Pressure,N*sizeof(double));
CopyToHost(Averages->Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages->Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages->Vel_z.get(),&Velocity[2*N],N*sizeof(double));
//...........................................................................
if (rank==0) printf("********************************************************\n");
@ -785,7 +785,7 @@ int main(int argc, char **argv)
timestep++;
// Run the analysis, blob identification, and write restart files
run_analysis(timestep,RESTART_INTERVAL,rank_info,Averages,last_ids,last_index,last_id_map,
run_analysis(timestep,RESTART_INTERVAL,rank_info,*Averages,last_ids,last_index,last_id_map,
Nx,Ny,Nz,pBC,beta,err,Phi,Pressure,Velocity,ID,f_even,f_odd,Den,
LocalRestartFile,tpool,work_ids);
@ -811,16 +811,16 @@ int main(int argc, char **argv)
// ************************************************************************
/* // Perform component averaging and write tcat averages
Averages.Initialize();
Averages.ComponentAverages();
Averages.SortBlobs();
Averages.PrintComponents(timestep);
Averages->Initialize();
Averages->ComponentAverages();
Averages->SortBlobs();
Averages->PrintComponents(timestep);
// ************************************************************************
int NumberComponents_NWP = ComputeGlobalPhaseComponent(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,Averages.PhaseID,1,Averages.Label_NWP);
int NumberComponents_NWP = ComputeGlobalPhaseComponent(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,Averages->PhaseID,1,Averages->Label_NWP);
printf("Number of non-wetting phase components: %i \n ",NumberComponents_NWP);
DeviceBarrier();
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages->Phase.get(),Phi,N*sizeof(double));
*/
// Create the MeshDataStruct
fillHalo<double> fillData(Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
@ -852,26 +852,26 @@ int main(int argc, char **argv)
BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(BlobIDVar);
fillData.copy(Averages.SDn,PhaseVar->data);
fillData.copy(Averages.SDs,SignDistVar->data);
fillData.copy(Averages.Label_NWP,BlobIDVar->data);
fillData.copy(Averages->SDn,PhaseVar->data);
fillData.copy(Averages->SDs,SignDistVar->data);
fillData.copy(Averages->Label_NWP,BlobIDVar->data);
IO::writeData( 0, meshData, 2 );
/* Averages.WriteSurfaces(0);
/* Averages->WriteSurfaces(0);
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
FILE *PHASE;
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Averages.Phase.get(),8,N,PHASE);
fwrite(Averages->Phase.get(),8,N,PHASE);
fclose(PHASE);
*/
/* sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
FILE *PRESS;
PRESS = fopen(LocalRankFilename,"wb");
fwrite(Averages.Press.get(),8,N,PRESS);
fwrite(Averages->Press.get(),8,N,PRESS);
fclose(PRESS);
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages->Phase.get(),Phi,N*sizeof(double));
double * Grad;
Grad = new double [3*N];
CopyToHost(Grad,ColorGrad,3*N*sizeof(double));