Modifying analysis in lbpm_color_simulator to run in seperate threads

This commit is contained in:
Mark Berrill
2015-08-22 11:03:46 -04:00
parent 4f6bceffc8
commit 3f771d9374
22 changed files with 5640 additions and 140 deletions

View File

@@ -27,6 +27,7 @@ extern void GlobalFlipInitD3Q19(double *dist_even, double *dist_odd, int Nx, int
X = Nx*nprocx;
Y = Ny*nprocy;
Z = Nz*nprocz;
NULL_USE(Z);
N = (Nx+2)*(Ny+2)*(Nz+2); // size of the array including halo
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
@@ -88,6 +89,7 @@ extern int GlobalCheckDebugDist(double *dist_even, double *dist_odd, int Nx, int
X = Nx*nprocx;
Y = Ny*nprocy;
Z = Nz*nprocz;
NULL_USE(Z);
N = (Nx+2)*(Ny+2)*(Nz+2); // size of the array including halo
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
@@ -161,7 +163,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
//**********************************
@@ -417,7 +418,6 @@ int main(int argc, char **argv)
starttime = MPI_Wtime();
//.........................................
sendtag = recvtag = 5;
//************ MAIN ITERATION LOOP (timing communications)***************************************/
while (timestep < 100){

View File

@@ -229,7 +229,7 @@ int main(int argc, char **argv)
// Variables that specify the computational domain
string FILENAME;
unsigned int nBlocks, nthreads;
//unsigned int nBlocks, nthreads;
int Nx,Ny,Nz; // local sub-domain size
int nspheres; // number of spheres in the packing
double Lx,Ly,Lz; // Domain length
@@ -248,8 +248,8 @@ int main(int argc, char **argv)
double fluid_isovalue,solid_isovalue;
fluid_isovalue = 0.0;
solid_isovalue = 0.0;
nBlocks = 32;
nthreads = 128;
//nBlocks = 32;
//nthreads = 128;
int RESTART_INTERVAL=20000;

View File

@@ -55,9 +55,6 @@ int main(int argc, char **argv)
int i,j,k,n;
// pmmc threshold values
double fluid_isovalue,solid_isovalue;
fluid_isovalue = 0.0;
solid_isovalue = 0.0;
if (rank==0){
//.......................................................................

View File

@@ -12,6 +12,9 @@
#include "common/MPI_Helpers.h"
#include "ProfilerApp.h"
#include "threadpool/thread_pool.h"
#include "lbpm_color_simulator.h"
//#define WRITE_SURFACES
@@ -108,7 +111,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
//**********************************
@@ -152,9 +154,9 @@ int main(int argc, char **argv)
int i,j,k;
// pmmc threshold values
double fluid_isovalue,solid_isovalue;
fluid_isovalue = 0.0;
solid_isovalue = 0.0;
//double fluid_isovalue,solid_isovalue;
//fluid_isovalue = 0.0;
//solid_isovalue = 0.0;
int RESTART_INTERVAL=20000;
@@ -255,8 +257,7 @@ int main(int argc, char **argv)
double Ps = -(das-dbs)/(das+dbs);
double rlxA = 1.f/tau;
double rlxB = 8.f*(2.f-rlxA)/(8.f-rlxA);
double xIntPos;
xIntPos = log((1.0+phi_s)/(1.0-phi_s))/(2.0*beta);
//double xIntPos = log((1.0+phi_s)/(1.0-phi_s))/(2.0*beta);
// Set the density values inside the solid based on the input value phi_s
das = (phi_s+1.0)*0.5;
@@ -299,6 +300,7 @@ int main(int argc, char **argv)
bool Restart;
if (InitialCondition==1) Restart=true;
else Restart=false;
NULL_USE(pBC); NULL_USE(velBC);
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
TwoPhase Averages(Dm);
@@ -523,11 +525,6 @@ int main(int argc, char **argv)
AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
//copies of data needed to perform checkpointing from cpu
double *cDen, *cDistEven, *cDistOdd;
cDen = new double[2*N];
cDistEven = new double[10*N];
cDistOdd = new double[9*N];
//...........................................................................
// Copy signed distance for device initialization
@@ -552,12 +549,18 @@ int main(int argc, char **argv)
if (Restart == true){
if (rank==0) printf("Reading restart file! \n");
// Read in the restart file to CPU buffers
double *cDen = new double[2*N];
double *cDistEven = new double[10*N];
double *cDistOdd = new double[9*N];
ReadCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
// Copy the restart data to the GPU
CopyToDevice(f_even,cDistEven,10*N*sizeof(double));
CopyToDevice(f_odd,cDistOdd,9*N*sizeof(double));
CopyToDevice(Den,cDen,2*N*sizeof(double));
DeviceBarrier();
delete [] cDen;
delete [] cDistEven;
delete [] cDistOdd;
MPI_Barrier(MPI_COMM_WORLD);
}
@@ -651,7 +654,6 @@ int main(int argc, char **argv)
CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
//...........................................................................
int timestep = -1;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
@@ -662,15 +664,21 @@ int main(int argc, char **argv)
starttime = MPI_Wtime();
//.........................................
sendtag = recvtag = 5;
// Copy the data to the CPU
err = 1.0;
double sat_w_previous = 1.01; // slightly impossible value!
if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
//************ MAIN ITERATION LOOP ***************************************/
PROFILE_START("Loop");
IntArray GlobalBlobID, GlobalBlobID2;
while (timestep < timestepMax && err > tol ){
int N_procs = ThreadPool::getNumberOfProcessors();
std::vector<int> procs(N_procs);
for (int i=0; i<N_procs; i++)
procs[i] = i;
ThreadPool::setProcessAffinity(procs);
int timestep = -1;
AnalysisWaitIdStruct work_ids;
ThreadPool tpool(2);
BlobIDstruct last_ids;
while (timestep < timestepMax && err > tol ) {
PROFILE_START("Update");
//*************************************************************************
@@ -773,98 +781,16 @@ int main(int argc, char **argv)
// Timestep completed!
timestep++;
//...................................................................
if (timestep%1000 == 995){
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Averages.Phase_tplus.get(),Phi,N*sizeof(double));
// Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
//...........................................................................
}
if (timestep%1000 == 0){
//...........................................................................
// Copy the data for for the analysis timestep
//...........................................................................
// Copy the phase from the GPU -> CPU
//...........................................................................
PROFILE_START("Copy phase");
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));
MPI_Barrier(MPI_COMM_WORLD);
PROFILE_STOP("Copy phase");
}
if ( timestep%1000 == 0){
// Compute the global blob id and compare to the previous version
PROFILE_START("Identify blobs and maps");
DeviceBarrier();
CopyToHost(Averages.Phase.get(),Phi,N*sizeof(double));
double vF = 0.0;
double vS = 0.0;
int nblobs2 = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,
Averages.Phase,Averages.SDs,vF,vS,GlobalBlobID2);
if ( !GlobalBlobID.empty() ) {
// Compute the timestep-timestep map
ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2);
// Renumber the current timestep's ids
}
Averages.NumberComponents_NWP = nblobs2;
Averages.Label_NWP.swap(GlobalBlobID2);
Averages.NumberComponents_WP = 1;
Averages.Label_WP.fill(0.0);
GlobalBlobID.swap(GlobalBlobID2);
PROFILE_STOP("Identify blobs and maps");
}
if (timestep%1000 == 5){
PROFILE_START("Compute dist");
//...........................................................................
// Copy the phase indicator field for the later timestep
DeviceBarrier();
CopyToHost(Averages.Phase_tminus.get(),Phi,N*sizeof(double));
// Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
//....................................................................
Averages.Initialize();
Averages.ComputeDelPhi();
Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
Averages.UpdateMeshValues();
Averages.ComputeLocal();
Averages.Reduce();
Averages.PrintAll(timestep);
Averages.Initialize();
Averages.ComponentAverages();
Averages.SortBlobs();
Averages.PrintComponents(timestep);
//....................................................................
PROFILE_STOP("Compute dist");
}
if (timestep%RESTART_INTERVAL == 0){
PROFILE_START("Save Checkpoint");
if (pBC){
//err = fabs(sat_w - sat_w_previous);
//sat_w_previous = sat_w;
if (rank==0) printf("Timestep %i: change in saturation since last checkpoint is %f \n", timestep, err);
}
else{
// Not clear yet
}
// Copy the data to the CPU
CopyToHost(cDistEven,f_even,10*N*sizeof(double));
CopyToHost(cDistOdd,f_odd,9*N*sizeof(double));
CopyToHost(cDen,Den,2*N*sizeof(double));
// Read in the restart file to CPU buffers
WriteCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
PROFILE_STOP("Save Checkpoint");
PROFILE_SAVE("lbpm_color_simulator",1);
}
// Run the analysis, blob identification, and write restart files
run_analysis(timestep,RESTART_INTERVAL,rank_info,Averages,last_ids,
Nx,Ny,Nz,pBC,beta,err,Phi,Pressure,Velocity,ID,f_even,f_odd,Den,
LocalRestartFile,tpool,work_ids);
}
tpool.wait_pool_finished();
PROFILE_STOP("Loop");
//************************************************************************/
//************************************************************************
DeviceBarrier();
MPI_Barrier(MPI_COMM_WORLD);
stoptime = MPI_Wtime();
@@ -881,14 +807,14 @@ int main(int argc, char **argv)
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
//************************************************************************/
// ************************************************************************
/* // Perform component averaging and write tcat averages
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);
printf("Number of non-wetting phase components: %i \n ",NumberComponents_NWP);
DeviceBarrier();
@@ -936,7 +862,7 @@ int main(int argc, char **argv)
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Averages.Phase.get(),8,N,PHASE);
fclose(PHASE);
*/
/* sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
FILE *PRESS;
PRESS = fopen(LocalRankFilename,"wb");
@@ -960,3 +886,5 @@ int main(int argc, char **argv)
MPI_Finalize();
// ****************************************************
}

View File

@@ -0,0 +1,220 @@
// Run the analysis, blob identification, and write restart files
enum AnalysisType{ AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02,
CopyAverages=0x02, CalcDist=0x02, CreateRestart=0x10 };
struct AnalysisWaitIdStruct {
ThreadPool::thread_id_t blobID;
ThreadPool::thread_id_t analysis;
ThreadPool::thread_id_t restart;
};
// Helper class to write the restart file from a seperate thread
class WriteRestartWorkItem: public ThreadPool::WorkItem
{
public:
WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cDen_,
std::shared_ptr<double> cDistEven_, std::shared_ptr<double>cDistOdd_, int N_ ):
filename(filename_), cDen(cDen_), cDistEven(cDistEven_), cDistOdd(cDistOdd_), N(N_) {}
virtual void run() {
PROFILE_START("Save Checkpoint",1);
WriteCheckpoint(filename,cDen.get(),cDistEven.get(),cDistOdd.get(),N);
PROFILE_STOP("Save Checkpoint",1);
PROFILE_SAVE("lbpm_color_simulator",1);
};
private:
WriteRestartWorkItem();
const char* filename;
std::shared_ptr<double> cDen, cDistEven, cDistOdd;
const int N;
};
// Helper class to compute the blob ids
typedef std::shared_ptr<std::pair<int,IntArray> > BlobIDstruct;
class BlobIdentificationWorkItem: public ThreadPool::WorkItem
{
public:
BlobIdentificationWorkItem( int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
std::shared_ptr<const double> phase_, const DoubleArray& dist_,
BlobIDstruct last_id_, BlobIDstruct new_id_ ):
Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), phase(phase_),
dist(dist_), last_id(last_id_), new_id(new_id_) { }
virtual void run() {
// Compute the global blob id and compare to the previous version
PROFILE_START("Identify blobs and maps",1);
double vF = 0.0;
double vS = 0.0;
new_id->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,
*phase,dist,vF,vS,new_id->second);
if ( last_id==NULL ) {
// Compute the timestep-timestep map
ID_map_struct map = computeIDMap(last_id->second,new_id->second);
// Renumber the current timestep's ids
}
PROFILE_STOP("Identify blobs and maps",1);
}
private:
BlobIdentificationWorkItem();
int Nx, Ny, Nz;
const RankInfoStruct& rank_info;
std::shared_ptr<const double> phase;
const DoubleArray& dist;
BlobIDstruct last_id, new_id;
};
// Helper class to run the analysis from within a thread
// Note: Averages will be modified after the constructor is called
class AnalysisWorkItem: public ThreadPool::WorkItem
{
public:
AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, BlobIDstruct ids, double beta_ ):
type(type_), timestep(timestep_), Averages(Averages_), blob_ids(ids), beta(beta_) { }
virtual void run() {
Averages.NumberComponents_NWP = blob_ids->first;
Averages.Label_NWP = blob_ids->second;
Averages.NumberComponents_WP = 1;
Averages.Label_WP.fill(0.0);
if ( (type&CopyPhaseIndicator) != 0 ) {
// Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
}
if ( (type&CalcDist) != 0 ) {
PROFILE_START("Compute dist",1);
// Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
Averages.Initialize();
Averages.ComputeDelPhi();
Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
Averages.UpdateMeshValues();
Averages.ComputeLocal();
Averages.Reduce();
Averages.PrintAll(timestep);
Averages.Initialize();
Averages.ComponentAverages();
Averages.SortBlobs();
Averages.PrintComponents(timestep);
PROFILE_STOP("Compute dist",1);
}
}
private:
AnalysisWorkItem();
AnalysisType type;
int timestep;
TwoPhase& Averages;
BlobIDstruct blob_ids;
double beta;
};
// Function to start the analysis
void run_analysis( int timestep, int restart_interval,
const RankInfoStruct& rank_info, TwoPhase& Averages, BlobIDstruct& last_ids,
int Nx, int Ny, int Nz, bool pBC, double beta, double err,
const double *Phi, double *Pressure, const double *Velocity,
const char *ID, const double *f_even, const double *f_odd, const double *Den,
const char *LocalRestartFile, ThreadPool& tpool, AnalysisWaitIdStruct& wait )
{
int N = Nx*Ny*Nz;
// Determin the analysis we want to perform
AnalysisType type = AnalyzeNone;
if ( timestep%1000 == 995 ) {
// Copy the phase indicator field for the earlier timestep
type = static_cast<AnalysisType>( type | CopyPhaseIndicator );
}
if ( timestep%1000 == 0 ) {
type = static_cast<AnalysisType>( type | CopyAverages );
type = static_cast<AnalysisType>( type | IdentifyBlobs );
}
if ( timestep%1000 == 5 ) {
type = static_cast<AnalysisType>( type | CalcDist );
}
if (timestep%restart_interval == 0) {
type = static_cast<AnalysisType>( type | CreateRestart );
}
// Return if we are not doing anything
if ( type == AnalyzeNone )
return;
PROFILE_START("start_analysis");
// Copy the appropriate variables to the host (so we can spawn new threads)
DeviceBarrier();
PROFILE_START("Copy data to host",1);
std::shared_ptr<double> phase;
if ( (type&CopyPhaseIndicator)!=0 || (type&CalcDist)!=0 ||
(type&CopyAverages)!=0 || (type&IdentifyBlobs)!=0 )
{
std::shared_ptr<double>(new double[N],DeleteArray<double>);
CopyToHost(phase.get(),Phi,N*sizeof(double));
}
if ( (type&CopyPhaseIndicator)!=0 ) {
memcpy(Averages.Phase_tplus.get(),phase.get(),N*sizeof(double));
}
if ( (type&CalcDist)!=0 ) {
memcpy(Averages.Phase_tminus.get(),phase.get(),N*sizeof(double));
}
if ( (type&CopyAverages) != 0 ) {
// Copy the members of Averages to the cpu (phase was copied above)
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
memcpy(Averages.Phase.get(),phase.get(),N*sizeof(double));
DeviceBarrier();
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));
}
std::shared_ptr<double> cDen, cDistEven, cDistOdd;
if ( (type&CreateRestart) != 0 ) {
// Copy restart data to the CPU
cDen = std::shared_ptr<double>(new double[2*N],DeleteArray<double>);
cDistEven = std::shared_ptr<double>(new double[10*N],DeleteArray<double>);
cDistOdd = std::shared_ptr<double>(new double[9*N],DeleteArray<double>);
CopyToHost(cDistEven.get(),f_even,10*N*sizeof(double));
CopyToHost(cDistOdd.get(),f_odd,9*N*sizeof(double));
CopyToHost(cDen.get(),Den,2*N*sizeof(double));
}
PROFILE_STOP("Copy data to host",1);
// Spawn threads to do blob identification work
if ( (type&IdentifyBlobs)!=0 ) {
BlobIDstruct new_ids;
ThreadPool::WorkItem *work = new BlobIdentificationWorkItem(
Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_ids);
work->add_dependency(wait.blobID);
last_ids = new_ids;
wait.blobID = tpool.add_work(work);
}
// Spawn threads to do the analysis work
if ( (type&CalcDist) != 0 ) {
ThreadPool::WorkItem *work = new AnalysisWorkItem(type,timestep,Averages,last_ids,beta);
work->add_dependency(wait.blobID);
work->add_dependency(wait.analysis);
wait.analysis = tpool.add_work(work);
}
// Spawn a thread to write the restart file
if ( (type&CreateRestart) != 0 ) {
int rank;
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
if (pBC) {
//err = fabs(sat_w - sat_w_previous);
//sat_w_previous = sat_w;
if (rank==0) printf("Timestep %i: change in saturation since last checkpoint is %f \n",timestep,err);
} else {
// Not clear yet
}
// Write the restart file (using a seperate thread)
WriteRestartWorkItem *work = new WriteRestartWorkItem(LocalRestartFile,cDen,cDistEven,cDistOdd,N);
work->add_dependency(wait.restart);
wait.restart = tpool.add_work(work);
}
PROFILE_STOP("start_analysis");
}

View File

@@ -31,7 +31,7 @@ inline void ReadDiscPacking(int ndiscs, double *List_cx, double *List_cy, double
char * line = new char[100];
// We will read until a blank like or end-of-file is reached
int count = 0;
while ( !feof(fid) && fgets(line,100,fid)>0 ) {
while ( !feof(fid) && fgets(line,100,fid)!=NULL ) {
char* line2 = line;
List_cx[count] = strtod(line2,&line2);
List_cy[count] = strtod(line2,&line2);
@@ -71,6 +71,7 @@ inline void SignedDistanceDiscPack(double *Distance, int ndiscs, double *List_cx
min_x = double(iproc*(Nx-2)-1)*hx;
min_y = double(jproc*(Ny-2)-1)*hy;
min_z = double(kproc*(Nz-2)-1)*hz;
NULL_USE(min_x);
//............................................
//............................................

View File

@@ -74,7 +74,6 @@ int main (int argc, char **argv)
DTMutableList<Point> local_nws_pts(20);
int n_local_nws_pts;
int n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg;
int c;
int newton_steps = 0;
//...........................................................................
@@ -153,10 +152,6 @@ int main (int argc, char **argv)
n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
n_nw_tris_beg = 0;// n_nw_tris;
n_ns_tris_beg = 0;//n_ns_tris;
n_ws_tris_beg = 0;//n_ws_tris;
// if there is a solid phase interface in the grid cell
if (Interface(SignDist,solid_isovalue,i,j,k) == 1){
@@ -387,9 +382,6 @@ int main (int argc, char **argv)
// Reset the triangle counts to zero
n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
n_nw_tris_beg = 0;// n_nw_tris;
// n_ns_tris_beg = 0;//n_ns_tris;
// n_ws_tris_beg = 0;//n_ws_tris;
// n_nws_seg_beg = n_nws_seg;
//*******************************************************************
}