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

Fixing bugs with thread based resturcturing of analysis in lbpm_color_simulator

Browse Source
This commit is contained in:
Mark Berrill 2015-08-24 14:11:26 -04:00
parent 6010fac583
commit 968e906032
3 changed files with 31 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
analysis/analysis.cpp
View File

@ -15,9 +15,7 @@ inline const TYPE* getPtr( const std::vector<TYPE>& x ) { return x.empty() ? NUL
int ComputeBlob( const Array<bool>& isPhase, IntArray& LocalBlobID, bool periodic, int start_id )
{
PROFILE_START("ComputeBlob",1);
ASSERT(isPhase.size(0)==LocalBlobID.size(0));
ASSERT(isPhase.size(1)==LocalBlobID.size(1));
ASSERT(isPhase.size(2)==LocalBlobID.size(2));
ASSERT(isPhase.size()==LocalBlobID.size());
const int Nx = isPhase.size(0); // maxima for the meshes
const int Ny = isPhase.size(1);
const int Nz = isPhase.size(2);
@ -136,9 +134,7 @@ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
double vF, double vS, IntArray& LocalBlobID, bool periodic )
{
PROFILE_START("ComputeLocalBlobIDs");
ASSERT(SignDist.size(0)==Phase.size(0));
ASSERT(SignDist.size(1)==Phase.size(1));
ASSERT(SignDist.size(2)==Phase.size(2));
ASSERT(SignDist.size()==Phase.size());
size_t Nx = Phase.size(0);
size_t Ny = Phase.size(1);
size_t Nz = Phase.size(2);
@ -170,7 +166,7 @@ int ComputeLocalPhaseComponent(const IntArray &PhaseID, int VALUE, IntArray &Com
size_t Nz = PhaseID.size(2);
size_t N = Nx*Ny*Nz;
// Compute the local blob ids
ComponentLabel.resize(Nx,Ny,Nz);
ComponentLabel.resize(Nx,Ny,Nz);
Array<bool> isPhase(Nx,Ny,Nz);
for (size_t i=0; i<N; i++) {
if ( PhaseID(i) == VALUE) {
@ -310,7 +306,7 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_
PROFILE_START("LocalToGlobalIDs",1);
const int rank = rank_info.rank[1][1][1];
int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
const int ngx = (IDs.size(0)-nx)/2;
const int ngy = (IDs.size(1)-ny)/2;
const int ngz = (IDs.size(2)-nz)/2;
@ -441,8 +437,6 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_inf
IntArray& GlobalBlobID )
{
PROFILE_START("ComputeGlobalBlobIDs");
int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
// First compute the local ids
int nblobs = ComputeLocalBlobIDs(Phase,SignDist,vF,vS,GlobalBlobID,false);
// Compute the global ids
@ -454,8 +448,6 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& r
const IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID )
{
PROFILE_START("ComputeGlobalPhaseComponent");
int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
// First compute the local ids
int nblobs = ComputeLocalPhaseComponent(PhaseID,VALUE,GlobalBlobID,false);
// Compute the global ids
@ -470,8 +462,8 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& r
******************************************************************/
void gatherSet( std::set<int>& set )
{
int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
std::vector<int> send_data(set.begin(),set.end());
int send_count = send_data.size();
std::vector<int> recv_count(nprocs,0), recv_disp(nprocs,0);
@ -487,8 +479,8 @@ void gatherSet( std::set<int>& set )
}
void gatherSrcIDMap( std::map<int,std::set<int> >& src_map )
{
int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
std::vector<int> send_data;
for (std::map<int,std::set<int> >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) {
int id = it->first;
@ -529,7 +521,7 @@ void addSrcDstIDs( int src_id, std::map<int,std::set<int> >& src_map,
}
ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 )
{
ASSERT(ID1.size(0)==ID2.size(0)&&ID1.size(1)==ID2.size(1)&&ID1.size(2)==ID2.size(2));
ASSERT(ID1.size()==ID2.size());
PROFILE_START("computeIDMap");
// Get a global list of all src/dst ids and the src map for each local blob

1
tests/lbpm_color_simulator.cpp
View File

@ -133,6 +133,7 @@ int main(int argc, char **argv)
//PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START("Main");
Utilities::setErrorHandlers();
// Variables that specify the computational domain
string FILENAME;

33
tests/lbpm_color_simulator.h
View File

@ -1,10 +1,12 @@
// Run the analysis, blob identification, and write restart files
#include "common/Array.h"
enum AnalysisType{ AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02,
CopyAverages=0x02, CalcDist=0x02, CreateRestart=0x10 };
// Structure used to store ids
struct AnalysisWaitIdStruct {
ThreadPool::thread_id_t blobID;
ThreadPool::thread_id_t analysis;
@ -20,10 +22,12 @@ public:
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() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress
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);
ThreadPool::WorkItem::d_state = 2; // Change state to finished
};
private:
WriteRestartWorkItem();
@ -39,30 +43,33 @@ 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_,
std::shared_ptr<const DoubleArray> 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() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress
// 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);
IntArray& ids = new_id->second;
new_id->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids);
if ( last_id==NULL ) {
// Compute the timestep-timestep map
ID_map_struct map = computeIDMap(last_id->second,new_id->second);
const IntArray& old_ids = new_id->second;
ID_map_struct map = computeIDMap(old_ids,ids);
// Renumber the current timestep's ids
}
PROFILE_STOP("Identify blobs and maps",1);
ThreadPool::WorkItem::d_state = 2; // Change state to finished
}
private:
BlobIdentificationWorkItem();
int Nx, Ny, Nz;
const RankInfoStruct& rank_info;
std::shared_ptr<const double> phase;
std::shared_ptr<const DoubleArray> phase;
const DoubleArray& dist;
BlobIDstruct last_id, new_id;
};
@ -76,6 +83,7 @@ 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() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress
Averages.NumberComponents_NWP = blob_ids->first;
Averages.Label_NWP = blob_ids->second;
Averages.NumberComponents_WP = 1;
@ -99,6 +107,7 @@ public:
Averages.PrintComponents(timestep);
PROFILE_STOP("Compute dist",1);
}
ThreadPool::WorkItem::d_state = 2; // Change state to finished
}
private:
AnalysisWorkItem();
@ -146,23 +155,23 @@ void run_analysis( int timestep, int restart_interval,
// 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;
std::shared_ptr<DoubleArray> 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));
phase = std::shared_ptr<DoubleArray>(new DoubleArray(Nx,Ny,Nz));
CopyToHost(phase->get(),Phi,N*sizeof(double));
}
if ( (type&CopyPhaseIndicator)!=0 ) {
memcpy(Averages.Phase_tplus.get(),phase.get(),N*sizeof(double));
memcpy(Averages.Phase_tplus.get(),phase->get(),N*sizeof(double));
}
if ( (type&CalcDist)!=0 ) {
memcpy(Averages.Phase_tminus.get(),phase.get(),N*sizeof(double));
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));
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));
@ -182,7 +191,7 @@ void run_analysis( int timestep, int restart_interval,
// Spawn threads to do blob identification work
if ( (type&IdentifyBlobs)!=0 ) {
BlobIDstruct new_ids;
BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray()));
ThreadPool::WorkItem *work = new BlobIdentificationWorkItem(
Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_ids);
work->add_dependency(wait.blobID);