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

Finishing id mapping between timesteps, modifying TwoPhase to print the new global id. Still testing.

Browse Source
This commit is contained in:
Mark Berrill 2015-08-25 14:32:36 -04:00
parent 968e906032
commit 63ef7d0fc6
15 changed files with 1792 additions and 1371 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
README.titan
View File

@ -24,8 +24,8 @@ cmake \
-D CMAKE_C_COMPILER:PATH=cc \ -D CMAKE_C_COMPILER:PATH=cc \
-D CMAKE_CXX_COMPILER:PATH=CC \ -D CMAKE_CXX_COMPILER:PATH=CC \
-D CMAKE_CXX_COMPILER:PATH=CC \ -D CMAKE_CXX_COMPILER:PATH=CC \
-D CMAKE_C_FLAGS="-DCBUB" \ -D CFLAGS="-DCBUB" \
-D CMAKE_CXX_FLAGS="-DCBUB" \ -D CXXFLAGS="-DCBUB" \
-D MPI_COMPILER:BOOL=TRUE \ -D MPI_COMPILER:BOOL=TRUE \
-D MPIEXEC=aprun \ -D MPIEXEC=aprun \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \ -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \

253
analysis/analysis.cpp
View File

@ -12,7 +12,7 @@ inline const TYPE* getPtr( const std::vector<TYPE>& x ) { return x.empty() ? NUL
/****************************************************************** /******************************************************************
* Compute the blobs * * Compute the blobs *
******************************************************************/ ******************************************************************/
int ComputeBlob( const Array<bool>& isPhase, IntArray& LocalBlobID, bool periodic, int start_id ) int ComputeBlob( const Array<bool>& isPhase, BlobIDArray& LocalBlobID, bool periodic, int start_id )
{ {
PROFILE_START("ComputeBlob",1); PROFILE_START("ComputeBlob",1);
ASSERT(isPhase.size()==LocalBlobID.size()); ASSERT(isPhase.size()==LocalBlobID.size());
@ -53,7 +53,7 @@ int ComputeBlob( const Array<bool>& isPhase, IntArray& LocalBlobID, bool periodi
std::vector<int> neighbor_ids; std::vector<int> neighbor_ids;
neighbor_ids.reserve(N_neighbors); neighbor_ids.reserve(N_neighbors);
const bool *isPhasePtr = isPhase.get(); const bool *isPhasePtr = isPhase.get();
int *LocalBlobIDPtr = LocalBlobID.get(); BlobIDType *LocalBlobIDPtr = LocalBlobID.get();
for (int z=0; z<Nz; z++) { for (int z=0; z<Nz; z++) {
for (int y=0; y<Ny; y++) { for (int y=0; y<Ny; y++) {
for (int x=0; x<Nx; x++) { for (int x=0; x<Nx; x++) {
@ -131,7 +131,7 @@ int ComputeBlob( const Array<bool>& isPhase, IntArray& LocalBlobID, bool periodi
* Compute the local blob ids * * Compute the local blob ids *
******************************************************************/ ******************************************************************/
int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist, int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
double vF, double vS, IntArray& LocalBlobID, bool periodic ) double vF, double vS, BlobIDArray& LocalBlobID, bool periodic )
{ {
PROFILE_START("ComputeLocalBlobIDs"); PROFILE_START("ComputeLocalBlobIDs");
ASSERT(SignDist.size()==Phase.size()); ASSERT(SignDist.size()==Phase.size());
@ -158,7 +158,7 @@ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
PROFILE_STOP("ComputeLocalBlobIDs"); PROFILE_STOP("ComputeLocalBlobIDs");
return nblobs; return nblobs;
} }
int ComputeLocalPhaseComponent(const IntArray &PhaseID, int VALUE, IntArray &ComponentLabel, bool periodic ) int ComputeLocalPhaseComponent(const IntArray &PhaseID, int VALUE, BlobIDArray &ComponentLabel, bool periodic )
{ {
PROFILE_START("ComputeLocalPhaseComponent"); PROFILE_START("ComputeLocalPhaseComponent");
size_t Nx = PhaseID.size(0); size_t Nx = PhaseID.size(0);
@ -186,7 +186,7 @@ int ComputeLocalPhaseComponent(const IntArray &PhaseID, int VALUE, IntArray &Com
/****************************************************************** /******************************************************************
* Reorder the global blob ids * * Reorder the global blob ids *
******************************************************************/ ******************************************************************/
static int ReorderBlobIDs2( IntArray& ID, int N_blobs, int ngx, int ngy, int ngz ) static int ReorderBlobIDs2( BlobIDArray& ID, int N_blobs, int ngx, int ngy, int ngz )
{ {
if ( N_blobs==0 ) if ( N_blobs==0 )
return 0; return 0;
@ -233,7 +233,7 @@ static int ReorderBlobIDs2( IntArray& ID, int N_blobs, int ngx, int ngy, int ngz
PROFILE_STOP("ReorderBlobIDs2",1); PROFILE_STOP("ReorderBlobIDs2",1);
return N_blobs2; return N_blobs2;
} }
void ReorderBlobIDs( IntArray& ID ) void ReorderBlobIDs( BlobIDArray& ID )
{ {
PROFILE_START("ReorderBlobIDs"); PROFILE_START("ReorderBlobIDs");
int tmp = ID.max()+1; int tmp = ID.max()+1;
@ -301,7 +301,7 @@ static bool updateLocalIds( const std::map<int64_t,int64_t>& remote_map,
return changed; return changed;
} }
static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info, static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info,
int nblobs, IntArray& IDs ) int nblobs, BlobIDArray& IDs )
{ {
PROFILE_START("LocalToGlobalIDs",1); PROFILE_START("LocalToGlobalIDs",1);
const int rank = rank_info.rank[1][1][1]; const int rank = rank_info.rank[1][1][1];
@ -327,9 +327,9 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_
if ( IDs(i) >= 0 ) if ( IDs(i) >= 0 )
IDs(i) += offset; IDs(i) += offset;
} }
const Array<int> LocalIDs = IDs; const BlobIDArray LocalIDs = IDs;
// Copy the ids and get the neighbors through the halos // Copy the ids and get the neighbors through the halos
fillHalo<int> fillData(rank_info,nx,ny,nz,1,1,1,0,1,true,true,true); fillHalo<BlobIDType> fillData(rank_info,nx,ny,nz,1,1,1,0,1,true,true,true);
fillData.fill(IDs); fillData.fill(IDs);
// Create a list of all neighbor ranks (excluding self) // Create a list of all neighbor ranks (excluding self)
std::vector<int> neighbors; std::vector<int> neighbors;
@ -414,14 +414,14 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_
for (size_t k=ngz; k<IDs.size(2)-ngz; k++) { for (size_t k=ngz; k<IDs.size(2)-ngz; k++) {
for (size_t j=ngy; j<IDs.size(1)-ngy; j++) { for (size_t j=ngy; j<IDs.size(1)-ngy; j++) {
for (size_t i=ngx; i<IDs.size(0)-ngx; i++) { for (size_t i=ngx; i<IDs.size(0)-ngx; i++) {
int id = IDs(i,j,k); BlobIDType id = IDs(i,j,k);
if ( id >= 0 ) if ( id >= 0 )
IDs(i,j,k) = final_map[id-offset]; IDs(i,j,k) = final_map[id-offset];
} }
} }
} }
// Fill the ghosts // Fill the ghosts
fillHalo<int> fillData2(rank_info,nx,ny,nz,1,1,1,0,1,true,true,true); fillHalo<BlobIDType> fillData2(rank_info,nx,ny,nz,1,1,1,0,1,true,true,true);
fillData2.fill(IDs); fillData2.fill(IDs);
// Reorder based on size (and compress the id space // Reorder based on size (and compress the id space
int N_blobs_global = ReorderBlobIDs2(IDs,N_blobs_tot,ngx,ngy,ngz); int N_blobs_global = ReorderBlobIDs2(IDs,N_blobs_tot,ngx,ngy,ngz);
@ -434,7 +434,7 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_
} }
int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info, int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info,
const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS, const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS,
IntArray& GlobalBlobID ) BlobIDArray& GlobalBlobID )
{ {
PROFILE_START("ComputeGlobalBlobIDs"); PROFILE_START("ComputeGlobalBlobIDs");
// First compute the local ids // First compute the local ids
@ -445,7 +445,7 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_inf
return nglobal; return nglobal;
} }
int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& rank_info, int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& rank_info,
const IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID ) const IntArray &PhaseID, int VALUE, BlobIDArray &GlobalBlobID )
{ {
PROFILE_START("ComputeGlobalPhaseComponent"); PROFILE_START("ComputeGlobalPhaseComponent");
// First compute the local ids // First compute the local ids
@ -460,34 +460,48 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& r
/****************************************************************** /******************************************************************
* Compute the mapping of blob ids between timesteps * * Compute the mapping of blob ids between timesteps *
******************************************************************/ ******************************************************************/
void gatherSet( std::set<int>& set ) template<class TYPE> inline MPI_Datatype getMPIType();
template<> inline MPI_Datatype getMPIType<int32_t>() { return MPI_INT; }
template<> inline MPI_Datatype getMPIType<int64_t>() {
if ( sizeof(int64_t)==sizeof(long int) )
return MPI_LONG;
else if ( sizeof(int64_t)==sizeof(double) )
return MPI_DOUBLE;
}
template<class TYPE>
void gatherSet( std::set<TYPE>& set )
{ {
int nprocs; int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs); MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
std::vector<int> send_data(set.begin(),set.end()); MPI_Datatype type = getMPIType<TYPE>();
std::vector<TYPE> send_data(set.begin(),set.end());
int send_count = send_data.size(); int send_count = send_data.size();
std::vector<int> recv_count(nprocs,0), recv_disp(nprocs,0); std::vector<int> recv_count(nprocs,0), recv_disp(nprocs,0);
MPI_Allgather(&send_count,1,MPI_INT,getPtr(recv_count),1,MPI_INT,MPI_COMM_WORLD); MPI_Allgather(&send_count,1,MPI_INT,getPtr(recv_count),1,MPI_INT,MPI_COMM_WORLD);
for (int i=1; i<nprocs; i++) for (int i=1; i<nprocs; i++)
recv_disp[i] = recv_disp[i-1] + recv_count[i-1]; recv_disp[i] = recv_disp[i-1] + recv_count[i-1];
std::vector<int> recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]); std::vector<TYPE> recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]);
MPI_Allgatherv(getPtr(send_data),send_count,MPI_INT, MPI_Allgatherv(getPtr(send_data),send_count,type,
getPtr(recv_data),getPtr(recv_count),getPtr(recv_disp),MPI_INT, getPtr(recv_data),getPtr(recv_count),getPtr(recv_disp),type,
MPI_COMM_WORLD); MPI_COMM_WORLD);
for (size_t i=0; i<recv_data.size(); i++) for (size_t i=0; i<recv_data.size(); i++)
set.insert(recv_data[i]); set.insert(recv_data[i]);
} }
void gatherSrcIDMap( std::map<int,std::set<int> >& src_map ) template<class TYPE>
void gatherSrcIDMap( std::map<TYPE,std::set<TYPE> >& src_map )
{ {
int nprocs; int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs); MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
std::vector<int> send_data; MPI_Datatype type = getMPIType<TYPE>();
for (std::map<int,std::set<int> >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) { std::vector<TYPE> send_data;
typename std::map<TYPE,std::set<TYPE> >::const_iterator it;
for (it=src_map.begin(); it!=src_map.end(); ++it) {
int id = it->first; int id = it->first;
const std::set<int>& src_ids = it->second; const std::set<TYPE>& src_ids = it->second;
send_data.push_back(id); send_data.push_back(id);
send_data.push_back(src_ids.size()); send_data.push_back(src_ids.size());
for (std::set<int>::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2) typename std::set<TYPE>::const_iterator it2;
for (it2=src_ids.begin(); it2!=src_ids.end(); ++it2)
send_data.push_back(*it2); send_data.push_back(*it2);
} }
int send_count = send_data.size(); int send_count = send_data.size();
@ -495,45 +509,45 @@ void gatherSrcIDMap( std::map<int,std::set<int> >& src_map )
MPI_Allgather(&send_count,1,MPI_INT,getPtr(recv_count),1,MPI_INT,MPI_COMM_WORLD); MPI_Allgather(&send_count,1,MPI_INT,getPtr(recv_count),1,MPI_INT,MPI_COMM_WORLD);
for (int i=1; i<nprocs; i++) for (int i=1; i<nprocs; i++)
recv_disp[i] = recv_disp[i-1] + recv_count[i-1]; recv_disp[i] = recv_disp[i-1] + recv_count[i-1];
std::vector<int> recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]); std::vector<TYPE> recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]);
MPI_Allgatherv(getPtr(send_data),send_count,MPI_INT, MPI_Allgatherv(getPtr(send_data),send_count,type,
getPtr(recv_data),getPtr(recv_count),getPtr(recv_disp),MPI_INT, getPtr(recv_data),getPtr(recv_count),getPtr(recv_disp),type,
MPI_COMM_WORLD); MPI_COMM_WORLD);
size_t i=0; size_t i=0;
while ( i < recv_data.size() ) { while ( i < recv_data.size() ) {
int id = recv_data[i]; int id = recv_data[i];
int count = recv_data[i+1]; int count = recv_data[i+1];
i += 2; i += 2;
std::set<int>& src_ids = src_map[id]; std::set<TYPE>& src_ids = src_map[id];
for (int j=0; j<count; j++,i++) for (int j=0; j<count; j++,i++)
src_ids.insert(recv_data[i]); src_ids.insert(recv_data[i]);
} }
} }
void addSrcDstIDs( int src_id, std::map<int,std::set<int> >& src_map, void addSrcDstIDs( BlobIDType src_id, std::map<BlobIDType,std::set<BlobIDType> >& src_map,
std::map<int,std::set<int> >& dst_map, std::set<int>& src, std::set<int>& dst ) std::map<BlobIDType,std::set<BlobIDType> >& dst_map, std::set<BlobIDType>& src, std::set<BlobIDType>& dst )
{ {
src.insert(src_id); src.insert(src_id);
const std::set<int>& dst_ids = dst_map[src_id]; const std::set<BlobIDType>& dst_ids = dst_map[src_id];
for (std::set<int>::const_iterator it=dst_ids.begin(); it!=dst_ids.end(); ++it) { for (std::set<BlobIDType>::const_iterator it=dst_ids.begin(); it!=dst_ids.end(); ++it) {
if ( dst.find(*it)==dst.end() ) if ( dst.find(*it)==dst.end() )
addSrcDstIDs(*it,dst_map,src_map,dst,src); addSrcDstIDs(*it,dst_map,src_map,dst,src);
} }
} }
ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 ) ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2 )
{ {
ASSERT(ID1.size()==ID2.size()); ASSERT(ID1.size()==ID2.size());
PROFILE_START("computeIDMap"); PROFILE_START("computeIDMap");
// Get a global list of all src/dst ids and the src map for each local blob // Get a global list of all src/dst ids and the src map for each local blob
std::set<int> src_set, dst_set; std::set<BlobIDType> src_set, dst_set;
std::map<int,std::set<int> > src_map; // Map of the src ids for each dst id std::map<BlobIDType,std::set<BlobIDType> > src_map; // Map of the src ids for each dst id
for (size_t i=0; i<ID1.length(); i++) { for (size_t i=0; i<ID1.length(); i++) {
if ( ID1(i)>=0 ) if ( ID1(i)>=0 )
src_set.insert(ID1(i)); src_set.insert(ID1(i));
if ( ID2(i)>=0 ) if ( ID2(i)>=0 )
dst_set.insert(ID2(i)); dst_set.insert(ID2(i));
if ( ID2(i)>=0 && ID1(i)>=0 ) { if ( ID2(i)>=0 && ID1(i)>=0 ) {
std::set<int>& src_ids = src_map[ID2(i)]; std::set<BlobIDType>& src_ids = src_map[ID2(i)];
src_ids.insert(ID1(i)); src_ids.insert(ID1(i));
} }
} }
@ -542,12 +556,12 @@ ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 )
gatherSet( dst_set ); gatherSet( dst_set );
gatherSrcIDMap( src_map ); gatherSrcIDMap( src_map );
// Compute the dst id map // Compute the dst id map
std::map<int,std::set<int> > dst_map; // Map of the dst ids for each src id std::map<BlobIDType,std::set<BlobIDType> > dst_map; // Map of the dst ids for each src id
for (std::map<int,std::set<int> >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) { for (std::map<BlobIDType,std::set<BlobIDType> >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) {
int id = it->first; BlobIDType id = it->first;
const std::set<int>& src_ids = it->second; const std::set<BlobIDType>& src_ids = it->second;
for (std::set<int>::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2) { for (std::set<BlobIDType>::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2) {
std::set<int>& dst_ids = dst_map[*it2]; std::set<BlobIDType>& dst_ids = dst_map[*it2];
dst_ids.insert(id); dst_ids.insert(id);
} }
} }
@ -555,53 +569,53 @@ ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 )
// Perform the mapping of ids // Perform the mapping of ids
ID_map_struct id_map; ID_map_struct id_map;
// Find new blobs // Find new blobs
for (std::set<int>::const_iterator it=dst_set.begin(); it!=dst_set.end(); ++it) { for (std::set<BlobIDType>::const_iterator it=dst_set.begin(); it!=dst_set.end(); ++it) {
if ( src_map.find(*it)==src_map.end() ) if ( src_map.find(*it)==src_map.end() )
id_map.created.push_back(*it); id_map.created.push_back(*it);
} }
// Fine blobs that disappeared // Fine blobs that disappeared
for (std::set<int>::const_iterator it=src_set.begin(); it!=src_set.end(); ++it) { for (std::set<BlobIDType>::const_iterator it=src_set.begin(); it!=src_set.end(); ++it) {
if ( dst_map.find(*it)==dst_map.end() ) if ( dst_map.find(*it)==dst_map.end() )
id_map.destroyed.push_back(*it); id_map.destroyed.push_back(*it);
} }
// Find blobs with a 1-to-1 mapping // Find blobs with a 1-to-1 mapping
std::vector<int> dst_list; std::vector<BlobIDType> dst_list;
dst_list.reserve(src_map.size()); dst_list.reserve(src_map.size());
for (std::map<int,std::set<int> >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) for (std::map<BlobIDType,std::set<BlobIDType> >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it)
dst_list.push_back(it->first); dst_list.push_back(it->first);
for (size_t i=0; i<dst_list.size(); i++) { for (size_t i=0; i<dst_list.size(); i++) {
int dst_id = dst_list[i]; int dst_id = dst_list[i];
const std::set<int>& src_ids = src_map[dst_id]; const std::set<BlobIDType>& src_ids = src_map[dst_id];
if ( src_ids.size()==1 ) { if ( src_ids.size()==1 ) {
int src_id = *src_ids.begin(); int src_id = *src_ids.begin();
const std::set<int>& dst_ids = dst_map[src_id]; const std::set<BlobIDType>& dst_ids = dst_map[src_id];
if ( dst_ids.size()==1 ) { if ( dst_ids.size()==1 ) {
ASSERT(*dst_ids.begin()==dst_id); ASSERT(*dst_ids.begin()==dst_id);
src_map.erase(dst_id); src_map.erase(dst_id);
dst_map.erase(src_id); dst_map.erase(src_id);
id_map.src_dst.push_back(std::pair<int,int>(src_id,dst_id)); id_map.src_dst.push_back(std::pair<BlobIDType,BlobIDType>(src_id,dst_id));
} }
} }
} }
// Handle merge/splits // Handle merge/splits
while ( !dst_map.empty() ) { while ( !dst_map.empty() ) {
// Get a lit of the src-dst ids // Get a lit of the src-dst ids
std::set<int> src, dst; std::set<BlobIDType> src, dst;
addSrcDstIDs( dst_map.begin()->first, src_map, dst_map, src, dst ); addSrcDstIDs( dst_map.begin()->first, src_map, dst_map, src, dst );
for (std::set<int>::const_iterator it=src.begin(); it!=src.end(); ++it) for (std::set<BlobIDType>::const_iterator it=src.begin(); it!=src.end(); ++it)
dst_map.erase(*it); dst_map.erase(*it);
for (std::set<int>::const_iterator it=dst.begin(); it!=dst.end(); ++it) for (std::set<BlobIDType>::const_iterator it=dst.begin(); it!=dst.end(); ++it)
src_map.erase(*it); src_map.erase(*it);
if ( src.size()==1 ) { if ( src.size()==1 ) {
// Bubble split // Bubble split
id_map.split.push_back( BlobIDSplitStruct(*src.begin(),std::vector<int>(dst.begin(),dst.end())) ); id_map.split.push_back( BlobIDSplitStruct(*src.begin(),std::vector<BlobIDType>(dst.begin(),dst.end())) );
} else if ( dst.size()==1 ) { } else if ( dst.size()==1 ) {
// Bubble merge // Bubble merge
id_map.merge.push_back( BlobIDMergeStruct(std::vector<int>(src.begin(),src.end()),*dst.begin()) ); id_map.merge.push_back( BlobIDMergeStruct(std::vector<BlobIDType>(src.begin(),src.end()),*dst.begin()) );
} else { } else {
// Bubble split/merge // Bubble split/merge
id_map.merge_split.push_back( BlobIDMergeSplitStruct( id_map.merge_split.push_back( BlobIDMergeSplitStruct(
std::vector<int>(src.begin(),src.end()), std::vector<int>(dst.begin(),dst.end() ) ) ); std::vector<BlobIDType>(src.begin(),src.end()), std::vector<BlobIDType>(dst.begin(),dst.end() ) ) );
} }
} }
ASSERT(src_map.empty()); ASSERT(src_map.empty());
@ -611,4 +625,135 @@ ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 )
} }
/******************************************************************
* Renumber the ids *
******************************************************************/
void getNewIDs( ID_map_struct& map, BlobIDType& id_max, std::vector<BlobIDType>& new_ids )
{
new_ids.resize(0);
// Renumber the ids that map directly
for (size_t i=0; i<map.src_dst.size(); i++) {
int id1 = map.src_dst[i].second;
int id2 = map.src_dst[i].first;
map.src_dst[i].second = id2;
if ( new_ids.size() < static_cast<size_t>(id1+1) )
new_ids.resize(id1+1,-1);
new_ids[id1] = id2;
}
// Renumber the created blobs to create new ids
for (size_t i=0; i<map.created.size(); i++) {
int id1 = map.created[i];
id_max++;
int id2 = id_max;
map.created[i] = id2;
if ( new_ids.size() < static_cast<size_t>(id1+1) )
new_ids.resize(id1+1,-1);
new_ids[id1] = id2;
}
// Renumber the blob splits to create new ids
for (size_t i=0; i<map.split.size(); i++) {
for (size_t j=0; j<map.split[i].second.size(); j++) {
int id1 = map.split[i].second[j];
id_max++;
int id2 = id_max;
map.split[i].second[j] = id2;
if ( new_ids.size() < static_cast<size_t>(id1+1) )
new_ids.resize(id1+1,-1);
new_ids[id1] = id2;
}
}
// Renumber the blob merges to create a new id
for (size_t i=0; i<map.merge.size(); i++) {
int id1 = map.merge[i].second;
id_max++;
int id2 = id_max;
map.merge[i].second = id2;
if ( new_ids.size() < static_cast<size_t>(id1+1) )
new_ids.resize(id1+1,-1);
new_ids[id1] = id2;
}
// Renumber the blob merge/splits to create new ids
for (size_t i=0; i<map.merge_split.size(); i++) {
for (size_t j=0; j<map.merge_split[i].second.size(); j++) {
int id1 = map.merge_split[i].second[j];
id_max++;
int id2 = id_max;
map.merge_split[i].second[j] = id2;
if ( new_ids.size() < static_cast<size_t>(id1+1) )
new_ids.resize(id1+1,-1);
new_ids[id1] = id2;
}
}
}
void renumberIDs( const std::vector<BlobIDType>& new_ids, BlobIDArray& IDs )
{
size_t N = IDs.length();
BlobIDType* ids = IDs.get();
for (size_t i=0; i<N; i++) {
BlobIDType id = ids[i];
if ( id>=0 )
ids[i] = new_ids[id];
}
}
/******************************************************************
* Write the id map for the given timestep *
******************************************************************/
void writeIDMap( const ID_map_struct& map, long long int timestep, const std::string& filename )
{
int rank;
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
if ( rank!=0 )
return;
bool empty = map.created.empty() && map.destroyed.empty() &&
map.split.empty() && map.merge.empty() && map.merge_split.empty();
for (size_t i=0; i<map.src_dst.size(); i++)
empty = empty && map.src_dst[i].first==map.src_dst[i].second;
if ( timestep!=0 && empty )
return;
FILE *fid = NULL;
if ( timestep==0 )
fid = fopen(filename.c_str(),"wb");
else
fid = fopen(filename.c_str(),"ab");
INSIST(fid!=NULL,std::string("Error opening file: ")+filename);
if ( empty ) {
fclose(fid);
return;
}
fprintf(fid,"%lli:",timestep);
for (size_t i=0; i<map.created.size(); i++)
fprintf(fid," -%lli",static_cast<long long int>(map.created[i]));
for (size_t i=0; i<map.destroyed.size(); i++)
fprintf(fid," %lli-",static_cast<long long int>(map.destroyed[i]));
for (size_t i=0; i<map.src_dst.size(); i++) {
if ( map.src_dst[i].first!=map.src_dst[i].second )
fprintf(fid," %lli-%lli",static_cast<long long int>(map.src_dst[i].first),
static_cast<long long int>(map.src_dst[i].second));
}
for (size_t i=0; i<map.split.size(); i++) {
fprintf(fid," %lli-%lli",static_cast<long long int>(map.split[i].first),
static_cast<long long int>(map.split[i].second[0]));
for (size_t j=1; j<map.split[i].second.size(); j++)
fprintf(fid,"/%lli",static_cast<long long int>(map.split[i].second[j]));
}
for (size_t i=0; i<map.merge.size(); i++) {
fprintf(fid," %lli",static_cast<long long int>(map.merge[i].first[0]));
for (size_t j=1; j<map.merge[i].first.size(); j++)
fprintf(fid,"/%lli",static_cast<long long int>(map.merge[i].first[j]));
fprintf(fid,"-%lli",static_cast<long long int>(map.merge[i].second));
}
for (size_t i=0; i<map.merge_split.size(); i++) {
fprintf(fid," %lli",static_cast<long long int>(map.merge_split[i].first[0]));
for (size_t j=1; j<map.merge_split[i].first.size(); j++)
fprintf(fid,"/%lli",static_cast<long long int>(map.merge_split[i].first[j]));
fprintf(fid,"-%lli",static_cast<long long int>(map.merge_split[i].second[0]));
for (size_t j=1; j<map.merge_split[i].second.size(); j++)
fprintf(fid,"/%lli",static_cast<long long int>(map.merge_split[i].second[j]));
}
fprintf(fid,"\n");
fclose(fid);
}

71
analysis/analysis.h
View File

@ -9,6 +9,11 @@
#include <vector> #include <vector>
// Define types to use for blob ids
typedef int32_t BlobIDType;
typedef Array<BlobIDType> BlobIDArray;
/*! /*!
* @brief Compute the blob * @brief Compute the blob
* @details Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob * @details Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
@ -22,7 +27,7 @@
* @return Returns the number of blobs * @return Returns the number of blobs
*/ */
int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist, int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
double vF, double vS, IntArray& LocalBlobID, bool periodic=true ); double vF, double vS, BlobIDArray& LocalBlobID, bool periodic=true );
/*! /*!
* @brief Compute blob of an arbitrary phase * @brief Compute blob of an arbitrary phase
@ -33,8 +38,7 @@ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
* @param[out] ComponentLabel * @param[out] ComponentLabel
* @param[in] periodic * @param[in] periodic
*/ */
int ComputeLocalPhaseComponent( const IntArray &PhaseID, int VALUE, IntArray &ComponentLabel, int ComputeLocalPhaseComponent( const IntArray &PhaseID, int VALUE, IntArray &ComponentLabel, bool periodic );
bool periodic );
/*! /*!
@ -54,7 +58,7 @@ int ComputeLocalPhaseComponent( const IntArray &PhaseID, int VALUE, IntArray &Co
*/ */
int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info, int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info,
const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS, const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS,
IntArray& GlobalBlobID ); BlobIDArray& GlobalBlobID );
/*! /*!
@ -71,7 +75,7 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_inf
* @return Return the number of components in the specified phase * @return Return the number of components in the specified phase
*/ */
int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& rank_info, int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& rank_info,
const IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID ); const IntArray &PhaseID, int VALUE, BlobIDArray &GlobalBlobID );
/*! /*!
@ -83,19 +87,26 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& r
* @param[in] nz Number of elements in the z-direction * @param[in] nz Number of elements in the z-direction
* @param[in/out] ID The ids of the blobs * @param[in/out] ID The ids of the blobs
*/ */
void ReorderBlobIDs( IntArray& ID ); void ReorderBlobIDs( BlobIDArray& ID );
typedef std::pair<int,std::vector<int> > BlobIDSplitStruct; typedef std::pair<BlobIDType,std::vector<BlobIDType> > BlobIDSplitStruct;
typedef std::pair<std::vector<int>,int> BlobIDMergeStruct; typedef std::pair<std::vector<BlobIDType>,BlobIDType> BlobIDMergeStruct;
typedef std::pair<std::vector<int>,std::vector<int> > BlobIDMergeSplitStruct; typedef std::pair<std::vector<BlobIDType>,std::vector<BlobIDType> > BlobIDMergeSplitStruct;
struct ID_map_struct { struct ID_map_struct {
std::vector<int> created; // list of new blobs that were created std::vector<BlobIDType> created; // list of new blobs that were created
std::vector<int> destroyed; // list of blobs that disappeared std::vector<BlobIDType> destroyed; // list of blobs that disappeared
std::vector<std::pair<int,int> > src_dst; // one-one mapping of blobs (first,second timestep id) std::vector<std::pair<BlobIDType,BlobIDType> > src_dst; // one-one mapping of blobs (first,second timestep id)
std::vector<BlobIDSplitStruct> split; // list of blobs that split std::vector<BlobIDSplitStruct> split; // list of blobs that split
std::vector<BlobIDMergeStruct> merge; // list of blobs that merged std::vector<BlobIDMergeStruct> merge; // list of blobs that merged
std::vector<BlobIDMergeSplitStruct> merge_split; // list of blobs that both merged and split std::vector<BlobIDMergeSplitStruct> merge_split; // list of blobs that both merged and split
//! Empty constructor
ID_map_struct() {}
//! Create initial map from N blobs (ordered 1:N-1)
ID_map_struct( int N ) {
created.resize(N);
for (int i=0; i<N; i++) { created[i]=i; }
}
}; };
@ -107,7 +118,41 @@ struct ID_map_struct {
* @param[in] ID1 The blob ids at the first timestep * @param[in] ID1 The blob ids at the first timestep
* @param[in] ID2 The blob ids at the second timestep * @param[in] ID2 The blob ids at the second timestep
*/ */
ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 ); ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2 );
/*!
* @brief Compute the new global ids based on the map
* @details This functions computes the time-consistent global ids for the
* current global id index
* @param[in/out] map The timestep mapping for the ids
* @param[in] id_max The globally largest id used previously
* @param[out] new_ids The newly renumbered blob ids (0:ids.max())
*/
void getNewIDs( ID_map_struct& map, BlobIDType& id_max, std::vector<BlobIDType>& new_ids );
/*!
* @brief Update the blob ids based on mapping
* @details This functions computes the map of blob ids between iterations.
* Note: we also update the map to reflect the new ids
* @param[out] new_ids The newly renumbered blob ids (0:ids.max())
* @param[in/out] IDs The blob ids to renumber
*/
void renumberIDs( const std::vector<BlobIDType>& new_id_list, BlobIDArray& IDs );
/*!
* @brief Write the ID map
* @details This functions writes the id map fo an iteration.
* If no ids changed, then nothing will be written
* Note: only rank 0 writes, and the file is created on timestep 0.
* @param[in] map The timestep mapping for the ids
* @param[in] timestep The current timestep (timestep 0 creates the file)
* @param[in] filename The filename to write/append
*/
void writeIDMap( const ID_map_struct& map, long long int timestep, const std::string& filename );
#endif #endif

11
cmake/libraries.cmake
View File

@ -221,10 +221,17 @@ MACRO( CONFIGURE_SYSTEM )
SET( SYSTEM_LDFLAGS ${SYSTEM_LDFLAGS} -rdynamic ) SET( SYSTEM_LDFLAGS ${SYSTEM_LDFLAGS} -rdynamic )
ENDIF() ENDIF()
ENDIF() ENDIF()
# Try to add -fPIC
SET( CMAKE_REQUIRED_FLAGS "${CMAKE_CXX_FLAGS} ${COVERAGE_FLAGS} -fPIC" )
CHECK_CXX_SOURCE_COMPILES( "int main() { return 0;}" fPIC )
IF ( fPIC )
SET( SYSTEM_LDFLAGS ${SYSTEM_LDFLAGS} -fPIC )
SET( SYSTEM_LDFLAGS ${SYSTEM_LDFLAGS} -fPIC )
ENDIF()
IF ( USING_GCC ) IF ( USING_GCC )
SET( SYSTEM_LIBS ${SYSTEM_LIBS} -lgfortran ) SET( SYSTEM_LIBS ${SYSTEM_LIBS} -lgfortran )
SET(CMAKE_C_FLAGS " ${CMAKE_C_FLAGS} -fPIC" ) SET(CFLAGS_EXTRA " ${CFLAGS_EXTRA} -fPIC" )
SET(CMAKE_CXX_FLAGS " ${CMAKE_CXX_FLAGS} -fPIC" ) SET(CXXFLAGS_EXTRA " ${CXXFLAGS_EXTRA} -fPIC" )
ENDIF() ENDIF()
ELSEIF( ${CMAKE_SYSTEM_NAME} STREQUAL "Darwin" ) ELSEIF( ${CMAKE_SYSTEM_NAME} STREQUAL "Darwin" )
# Max specific system libraries # Max specific system libraries

10
common/Domain.h
View File

@ -16,7 +16,7 @@
#include "common/MPI_Helpers.h" #include "common/MPI_Helpers.h"
#include "common/Communication.h" #include "common/Communication.h"
int MAX_BLOB_COUNT=50; static int MAX_BLOB_COUNT=50;
using namespace std; using namespace std;
@ -211,7 +211,7 @@ static inline void fgetl( char * str, int num, FILE * stream )
if ( 0 ) {char *temp = (char *)&ptr; temp++;} if ( 0 ) {char *temp = (char *)&ptr; temp++;}
} }
Domain::~Domain(){ inline Domain::~Domain(){
delete sendData_x; delete sendData_x;
delete sendData_y; delete sendData_y;
delete sendData_z; delete sendData_z;
@ -250,7 +250,7 @@ Domain::~Domain(){
delete recvData_YZ; delete recvData_YZ;
} }
void Domain::InitializeRanks() inline void Domain::InitializeRanks()
{ {
// map the rank to the block index // map the rank to the block index
iproc = rank%nprocx; iproc = rank%nprocx;
@ -282,7 +282,7 @@ void Domain::InitializeRanks()
} }
void Domain::CommInit(MPI_Comm Communicator){ inline void Domain::CommInit(MPI_Comm Communicator){
int i,j,k,n; int i,j,k,n;
int sendtag = 21; int sendtag = 21;
int recvtag = 21; int recvtag = 21;
@ -696,7 +696,7 @@ inline void Domain::CommunicateMeshHalo(DoubleArray &Mesh)
UnpackMeshData(recvList_YZ, recvCount_YZ ,recvData_YZ, MeshData); UnpackMeshData(recvList_YZ, recvCount_YZ ,recvData_YZ, MeshData);
} }
void Domain::BlobComm(MPI_Comm Communicator){ inline void Domain::BlobComm(MPI_Comm Communicator){
//...................................................................................... //......................................................................................
int sendtag, recvtag; int sendtag, recvtag;
sendtag = recvtag = 51; sendtag = recvtag = 51;

1392
common/TwoPhase.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

1262
common/TwoPhase.h
View File

File diff suppressed because it is too large Load Diff

16
common/pmmc.h
View File

@ -1,3 +1,7 @@
#ifndef pmmc_INC
#define pmmc_INC
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <iostream> #include <iostream>
@ -10,7 +14,7 @@
using namespace std; using namespace std;
int edgeTable[256]={ static int edgeTable[256]={
0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
@ -43,8 +47,8 @@ int edgeTable[256]={
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 }; 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 };
char triTable[256][16] = static char triTable[256][16] =
{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
@ -735,7 +739,7 @@ inline bool Solid( DoubleArray &A, int i, int j, int k){
return X; return X;
} }
//------------------------------------------------------------------------------- //-------------------------------------------------------------------------------
Point VertexInterp(const Point &p1, const Point &p2, double valp1, double valp2) inline Point VertexInterp(const Point &p1, const Point &p2, double valp1, double valp2)
{ {
return (p1 + (-valp1 / (valp2 - valp1)) * (p2 - p1)); return (p1 + (-valp1 / (valp2 - valp1)) * (p2 - p1));
} }
@ -4032,3 +4036,7 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
//............................................................................. //.............................................................................
} }
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
#endif

26
pmmc/pmmc.cpp
View File

@ -11,7 +11,7 @@
using namespace std; using namespace std;
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
inline int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator, int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator,
DoubleArray &F, DoubleArray &S, double vf, double vs, int startx, int starty, DoubleArray &F, DoubleArray &S, double vf, double vs, int startx, int starty,
int startz, IntArray &temp) int startz, IntArray &temp)
{ {
@ -140,7 +140,7 @@ inline int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indi
} }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
/* /*
inline DoubleArray SOLVE( DoubleArray &A, DoubleArray &b) DoubleArray SOLVE( DoubleArray &A, DoubleArray &b)
{ {
// solves the system A*x = b exactly // solves the system A*x = b exactly
@ -158,7 +158,7 @@ inline DoubleArray SOLVE( DoubleArray &A, DoubleArray &b)
return solution; return solution;
} }
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
inline Point NEWTON(Point x0, DoubleArray &F, double &v,DoubleArray &S, int i, int j, int k, int &newton_steps) Point NEWTON(Point x0, DoubleArray &F, double &v,DoubleArray &S, int i, int j, int k, int &newton_steps)
{ {
Point pt; Point pt;
// Performs a Newton iteration to compute the point x from initial guess x0 // Performs a Newton iteration to compute the point x from initial guess x0
@ -321,7 +321,7 @@ inline Point NEWTON(Point x0, DoubleArray &F, double &v,DoubleArray &S, int i,
*/ */
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
inline bool vertexcheck(Point &P, int n, int pos, DTMutableList<Point> &cellvertices){ bool vertexcheck(Point &P, int n, int pos, DTMutableList<Point> &cellvertices){
// returns true if P is a new vertex (one previously unencountered // returns true if P is a new vertex (one previously unencountered
bool V = 1; bool V = 1;
@ -337,7 +337,7 @@ inline bool vertexcheck(Point &P, int n, int pos, DTMutableList<Point> &cellvert
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
inline bool ShareSide( Point &A, Point &B) bool ShareSide( Point &A, Point &B)
{ {
// returns true if points A and B share an x,y, or z coordinate // returns true if points A and B share an x,y, or z coordinate
bool l = 0; bool l = 0;
@ -360,7 +360,7 @@ inline bool ShareSide( Point &A, Point &B)
} }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
inline bool Interface( DoubleArray &A, const double v, int i, int j, int k){ bool Interface( DoubleArray &A, const double v, int i, int j, int k){
// returns true if grid cell i, j, k contains a section of the interface // returns true if grid cell i, j, k contains a section of the interface
bool Y = 0; bool Y = 0;
@ -416,7 +416,7 @@ inline bool Interface( DoubleArray &A, const double v, int i, int j, int k){
} }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
inline bool Fluid_Interface( DoubleArray &A, DoubleArray &S, const double v, int i, int j, int k){ bool Fluid_Interface( DoubleArray &A, DoubleArray &S, const double v, int i, int j, int k){
// returns true if grid cell i, j, k contains a section of the interface // returns true if grid cell i, j, k contains a section of the interface
bool Y = 0; bool Y = 0;
@ -470,7 +470,7 @@ inline bool Fluid_Interface( DoubleArray &A, DoubleArray &S, const double v, in
return Y; return Y;
} }
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
inline bool Solid( DoubleArray &A, int i, int j, int k){ bool Solid( DoubleArray &A, int i, int j, int k){
bool X = 0; bool X = 0;
@ -504,7 +504,7 @@ inline bool Solid( DoubleArray &A, int i, int j, int k){
} }
//------------------------------------------------------------------------------- //-------------------------------------------------------------------------------
//------------------------------------------------------------------------------- //-------------------------------------------------------------------------------
inline void SOL_SURF(DoubleArray &A, const double &v, DoubleArray &B, const double &isovalue, void SOL_SURF(DoubleArray &A, const double &v, DoubleArray &B, const double &isovalue,
int i,int j,int k, int m, int n, int o, DTMutableList<Point> int i,int j,int k, int m, int n, int o, DTMutableList<Point>
&cellvertices, int &lengthvertices, IntArray &Tlist, int &nTris, &cellvertices, int &lengthvertices, IntArray &Tlist, int &nTris,
DoubleArray &values){ DoubleArray &values){
@ -941,7 +941,7 @@ inline void SOL_SURF(DoubleArray &A, const double &v, DoubleArray &B, const doub
} }
} }
//------------------------------------------------------------------------------- //-------------------------------------------------------------------------------
inline void TRIM(DTMutableList<Point> &local_sol_pts, int &n_local_sol_pts, double isovalue, void TRIM(DTMutableList<Point> &local_sol_pts, int &n_local_sol_pts, double isovalue,
IntArray &local_sol_tris, int &n_local_sol_tris, IntArray &local_sol_tris, int &n_local_sol_tris,
DTMutableList<Point> &ns_pts, int &n_ns_pts, IntArray &ns_tris, DTMutableList<Point> &ns_pts, int &n_ns_pts, IntArray &ns_tris,
int &n_ns_tris, DTMutableList<Point> &ws_pts, int &n_ws_pts, int &n_ns_tris, DTMutableList<Point> &ws_pts, int &n_ws_pts,
@ -1374,7 +1374,7 @@ inline void TRIM(DTMutableList<Point> &local_sol_pts, int &n_local_sol_pts, doub
} }
} }
//------------------------------------------------------------------------------- //-------------------------------------------------------------------------------
inline void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k, void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k,
DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris, DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris,
int &n_nw_tris) int &n_nw_tris)
{ {
@ -2128,7 +2128,7 @@ else if ( A(i,j+1,k+1) == 0){
} }
//------------------------------------------------------------------------------- //-------------------------------------------------------------------------------
inline void EDGE(DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k, int &m, int &n, int &o, void EDGE(DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k, int &m, int &n, int &o,
DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris, int &n_nw_tris, DTMutableList<Point> &nw_pts, int &n_nw_pts, IntArray &nw_tris, int &n_nw_tris,
DTMutableList<Point> &local_nws_pts, int &n_local_nws_pts) DTMutableList<Point> &local_nws_pts, int &n_local_nws_pts)
{ {
@ -2491,7 +2491,7 @@ inline void EDGE(DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j,
} }
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
inline void ComputeAreasPMMC(IntArray &cubeList, int start, int finish, void ComputeAreasPMMC(IntArray &cubeList, int start, int finish,
DoubleArray &F, DoubleArray &S, double vF, double vS, DoubleArray &F, DoubleArray &S, double vF, double vS,
double &blob_volume, double &ans, double &aws, double &awn, double &lwns, double &blob_volume, double &ans, double &aws, double &awn, double &lwns,
int Nx, int Ny, int Nz) int Nx, int Ny, int Nz)

2
tests/CMakeLists.txt
View File

@ -27,9 +27,9 @@ ADD_LBPM_TEST( TestInterfaceSpeed )
ADD_LBPM_TEST( TestSphereCurvature ) ADD_LBPM_TEST( TestSphereCurvature )
ADD_LBPM_TEST( TestContactAngle ) ADD_LBPM_TEST( TestContactAngle )
ADD_LBPM_TEST_1_2_4( TestTwoPhase ) ADD_LBPM_TEST_1_2_4( TestTwoPhase )
ADD_LBPM_TEST_1_2_4( TestBlobIdentify )
ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 ) ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 )
ADD_LBPM_TEST_PARALLEL( TestBlobAnalyze 8 ) ADD_LBPM_TEST_PARALLEL( TestBlobAnalyze 8 )
ADD_LBPM_TEST_PARALLEL( TestBlobIdentify 1 )
ADD_LBPM_TEST_PARALLEL( TestSegDist 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_PARALLEL( TestMassConservationD3Q7 1 )

32
tests/TestBlobIdentify.cpp
View File

@ -22,6 +22,16 @@ inline double rand2()
} }
// Test if all ranks agree on a value
bool allAgree( int x ) {
int min, max;
MPI_Allreduce(&x,&min,1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
MPI_Allreduce(&x,&max,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
return min==max;
}
// Structure to hold a bubble
struct bubble_struct { struct bubble_struct {
Point center; Point center;
double radius; double radius;
@ -209,6 +219,7 @@ int main(int argc, char **argv)
fillData.copy(SignDist,SignDistVar->data); fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data); fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( 0, meshData, 2 ); IO::writeData( 0, meshData, 2 );
writeIDMap(ID_map_struct(nblobs),0,"lbpm_id_map.txt");
int save_it = 1; int save_it = 1;
// Check the results // Check the results
@ -222,6 +233,7 @@ int main(int argc, char **argv)
// Move the blobs and connect them to the previous ids // Move the blobs and connect them to the previous ids
PROFILE_START("constant velocity test"); PROFILE_START("constant velocity test");
if ( rank==0 ) { printf("Running constant velocity blob test\n"); } if ( rank==0 ) { printf("Running constant velocity blob test\n"); }
int id_max = nblobs-1;
for (int i=0; i<20; i++, save_it++) { for (int i=0; i<20; i++, save_it++) {
// Shift all the data // Shift all the data
shift_data( Phase, 3, -2, 1, rank_info ); shift_data( Phase, 3, -2, 1, rank_info );
@ -236,6 +248,11 @@ int main(int argc, char **argv)
} }
// Identify the blob maps and renumber the ids // Identify the blob maps and renumber the ids
ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2); ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2);
std::swap(GlobalBlobID,GlobalBlobID2);
std::vector<BlobIDType> new_list;
getNewIDs(map,id_max,new_list);
renumberIDs(new_list,GlobalBlobID);
writeIDMap(map,save_it,"lbpm_id_map.txt");
bool pass = (int)map.src_dst.size()==nblobs; bool pass = (int)map.src_dst.size()==nblobs;
pass = pass && map.created.empty(); pass = pass && map.created.empty();
pass = pass && map.destroyed.empty(); pass = pass && map.destroyed.empty();
@ -244,11 +261,11 @@ int main(int argc, char **argv)
pass = pass && map.split.empty(); pass = pass && map.split.empty();
pass = pass && map.merge.empty(); pass = pass && map.merge.empty();
pass = pass && map.merge_split.empty(); pass = pass && map.merge_split.empty();
pass = pass && id_max==nblobs-1;
if ( !pass ) { if ( !pass ) {
printf("Error, blob ids do not match in constant velocity test\n"); printf("Error, blob ids do not match in constant velocity test\n");
N_errors++; N_errors++;
} }
GlobalBlobID = GlobalBlobID2;
// Save the results // Save the results
fillData.copy(Phase,PhaseVar->data); fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data); fillData.copy(SignDist,SignDistVar->data);
@ -279,6 +296,7 @@ int main(int argc, char **argv)
fillData.copy(GlobalBlobID,BlobIDVar->data); fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( save_it, meshData, 2 ); IO::writeData( save_it, meshData, 2 );
save_it++; save_it++;
id_max = nblobs-1;
for (int i=0; i<25; i++, save_it++) { for (int i=0; i<25; i++, save_it++) {
// Move the bubbles // Move the bubbles
for (size_t j=0; j<bubbles.size(); j++) { for (size_t j=0; j<bubbles.size(); j++) {
@ -294,6 +312,17 @@ int main(int argc, char **argv)
int nblobs2 = ComputeGlobalBlobIDs(nx,ny,nz,rank_info,Phase,SignDist,vF,vS,GlobalBlobID2); int nblobs2 = ComputeGlobalBlobIDs(nx,ny,nz,rank_info,Phase,SignDist,vF,vS,GlobalBlobID2);
// Identify the blob maps and renumber the ids // Identify the blob maps and renumber the ids
ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2); ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2);
std::swap(GlobalBlobID,GlobalBlobID2);
std::vector<BlobIDType> new_list;
getNewIDs(map,id_max,new_list);
renumberIDs(new_list,GlobalBlobID);
writeIDMap(map,save_it,"lbpm_id_map.txt");
if ( !allAgree(id_max) ) {
if ( rank==0 )
printf("All ranks do not agree on id_max\n");
N_errors++;
break;
}
int N1 = 0; int N1 = 0;
int N2 = 0; int N2 = 0;
if ( rank==0 ) { if ( rank==0 ) {
@ -350,7 +379,6 @@ int main(int argc, char **argv)
N_errors++; N_errors++;
} }
nblobs = nblobs2; nblobs = nblobs2;
GlobalBlobID = GlobalBlobID2;
// Save the results // Save the results
fillData.copy(Phase,PhaseVar->data); fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data); fillData.copy(SignDist,SignDistVar->data);

8
tests/lbpm_color_simulator.cpp
View File

@ -677,8 +677,10 @@ int main(int argc, char **argv)
ThreadPool::setProcessAffinity(procs); ThreadPool::setProcessAffinity(procs);
int timestep = -1; int timestep = -1;
AnalysisWaitIdStruct work_ids; AnalysisWaitIdStruct work_ids;
ThreadPool tpool(2); ThreadPool tpool(0);
BlobIDstruct last_ids; BlobIDstruct last_ids, last_index;
BlobIDList last_id_map;
writeIDMap(ID_map_struct(),0,id_map_filename);
while (timestep < timestepMax && err > tol ) { while (timestep < timestepMax && err > tol ) {
PROFILE_START("Update"); PROFILE_START("Update");
@ -784,7 +786,7 @@ int main(int argc, char **argv)
timestep++; timestep++;
// Run the analysis, blob identification, and write restart files // Run the analysis, blob identification, and write restart files
run_analysis(timestep,RESTART_INTERVAL,rank_info,Averages,last_ids, 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, Nx,Ny,Nz,pBC,beta,err,Phi,Pressure,Velocity,ID,f_even,f_odd,Den,
LocalRestartFile,tpool,work_ids); LocalRestartFile,tpool,work_ids);

56
tests/lbpm_color_simulator.h
View File

@ -38,40 +38,54 @@ private:
// Helper class to compute the blob ids // Helper class to compute the blob ids
static const std::string id_map_filename = "lbpm_id_map.txt";
typedef std::shared_ptr<std::pair<int,IntArray> > BlobIDstruct; typedef std::shared_ptr<std::pair<int,IntArray> > BlobIDstruct;
typedef std::shared_ptr<std::vector<BlobIDType> > BlobIDList;
class BlobIdentificationWorkItem: public ThreadPool::WorkItem class BlobIdentificationWorkItem: public ThreadPool::WorkItem
{ {
public: public:
BlobIdentificationWorkItem( int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, BlobIdentificationWorkItem( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_, std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
BlobIDstruct last_id_, BlobIDstruct new_id_ ): BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ):
Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), phase(phase_), timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
dist(dist_), last_id(last_id_), new_id(new_id_) { } phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_) { }
virtual void run() { virtual void run() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress ThreadPool::WorkItem::d_state = 1; // Change state to in progress
// Compute the global blob id and compare to the previous version // Compute the global blob id and compare to the previous version
PROFILE_START("Identify blobs and maps",1); PROFILE_START("Identify blobs and maps",1);
double vF = 0.0; double vF = 0.0;
double vS = 0.0; double vS = 0.0;
IntArray& ids = new_id->second; IntArray& ids = new_index->second;
new_id->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids); new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids);
if ( last_id==NULL ) { static int max_id = -1;
new_id->first = new_index->first;
new_id->second = new_index->second;
if ( last_id!=NULL ) {
// Compute the timestep-timestep map // Compute the timestep-timestep map
const IntArray& old_ids = new_id->second; const IntArray& old_ids = last_id->second;
ID_map_struct map = computeIDMap(old_ids,ids); ID_map_struct map = computeIDMap(old_ids,ids);
// Renumber the current timestep's ids // Renumber the current timestep's ids
getNewIDs(map,max_id,*new_list);
renumberIDs(*new_list,new_id->second);
writeIDMap(map,timestep,id_map_filename);
} else {
max_id = -1;
ID_map_struct map(new_id->first);
getNewIDs(map,max_id,*new_list);
writeIDMap(map,timestep,id_map_filename);
} }
PROFILE_STOP("Identify blobs and maps",1); PROFILE_STOP("Identify blobs and maps",1);
ThreadPool::WorkItem::d_state = 2; // Change state to finished ThreadPool::WorkItem::d_state = 2; // Change state to finished
} }
private: private:
BlobIdentificationWorkItem(); BlobIdentificationWorkItem();
int timestep;
int Nx, Ny, Nz; int Nx, Ny, Nz;
const RankInfoStruct& rank_info; const RankInfoStruct& rank_info;
std::shared_ptr<const DoubleArray> phase; std::shared_ptr<const DoubleArray> phase;
const DoubleArray& dist; const DoubleArray& dist;
BlobIDstruct last_id, new_id; BlobIDstruct last_id, new_index, new_id;
BlobIDList new_list;
}; };
@ -80,12 +94,15 @@ private:
class AnalysisWorkItem: public ThreadPool::WorkItem class AnalysisWorkItem: public ThreadPool::WorkItem
{ {
public: public:
AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, BlobIDstruct ids, double beta_ ): AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_,
type(type_), timestep(timestep_), Averages(Averages_), blob_ids(ids), beta(beta_) { } BlobIDstruct ids, BlobIDList id_list_, double beta_ ):
type(type_), timestep(timestep_), Averages(Averages_),
blob_ids(ids), id_list(id_list_), beta(beta_) { }
virtual void run() { virtual void run() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress ThreadPool::WorkItem::d_state = 1; // Change state to in progress
Averages.NumberComponents_NWP = blob_ids->first; Averages.NumberComponents_NWP = blob_ids->first;
Averages.Label_NWP = blob_ids->second; Averages.Label_NWP = blob_ids->second;
Averages.Label_NWP_map = *id_list;
Averages.NumberComponents_WP = 1; Averages.NumberComponents_WP = 1;
Averages.Label_WP.fill(0.0); Averages.Label_WP.fill(0.0);
if ( (type&CopyPhaseIndicator) != 0 ) { if ( (type&CopyPhaseIndicator) != 0 ) {
@ -115,13 +132,15 @@ private:
int timestep; int timestep;
TwoPhase& Averages; TwoPhase& Averages;
BlobIDstruct blob_ids; BlobIDstruct blob_ids;
BlobIDList id_list;
double beta; double beta;
}; };
// Function to start the analysis // Function to start the analysis
void run_analysis( int timestep, int restart_interval, void run_analysis( int timestep, int restart_interval,
const RankInfoStruct& rank_info, TwoPhase& Averages, BlobIDstruct& last_ids, const RankInfoStruct& rank_info, TwoPhase& Averages,
BlobIDstruct& last_ids, BlobIDstruct& last_index, BlobIDList& last_id_map,
int Nx, int Ny, int Nz, bool pBC, double beta, double err, int Nx, int Ny, int Nz, bool pBC, double beta, double err,
const double *Phi, double *Pressure, const double *Velocity, const double *Phi, double *Pressure, const double *Velocity,
const char *ID, const double *f_even, const double *f_odd, const double *Den, const char *ID, const double *f_even, const double *f_odd, const double *Den,
@ -191,17 +210,22 @@ void run_analysis( int timestep, int restart_interval,
// Spawn threads to do blob identification work // Spawn threads to do blob identification work
if ( (type&IdentifyBlobs)!=0 ) { if ( (type&IdentifyBlobs)!=0 ) {
BlobIDstruct new_index(new std::pair<int,IntArray>(0,IntArray()));
BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray())); BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray()));
ThreadPool::WorkItem *work = new BlobIdentificationWorkItem( BlobIDList new_list(new std::vector<BlobIDType>());
Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_ids); ThreadPool::WorkItem *work = new BlobIdentificationWorkItem(timestep,
Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_index,new_ids,new_list);
work->add_dependency(wait.blobID); work->add_dependency(wait.blobID);
last_index = new_index;
last_ids = new_ids; last_ids = new_ids;
last_id_map = new_list;
wait.blobID = tpool.add_work(work); wait.blobID = tpool.add_work(work);
} }
// Spawn threads to do the analysis work // Spawn threads to do the analysis work
if ( (type&CalcDist) != 0 ) { if ( (type&CalcDist) != 0 ) {
ThreadPool::WorkItem *work = new AnalysisWorkItem(type,timestep,Averages,last_ids,beta); ThreadPool::WorkItem *work = new AnalysisWorkItem(
type,timestep,Averages,last_index,last_id_map,beta);
work->add_dependency(wait.blobID); work->add_dependency(wait.blobID);
work->add_dependency(wait.analysis); work->add_dependency(wait.analysis);
wait.analysis = tpool.add_work(work); wait.analysis = tpool.add_work(work);

18
visit/avtLBMFileFormat.C
View File

@ -99,9 +99,23 @@
#endif #endif
// Null stream buffer
class NullBufferClass : public std::streambuf {
public:
virtual std::streamsize xsputn(const char *s, std::streamsize n) { return n; }
virtual int overflow(int c) { return 1; }
virtual void setOutputStream( std::ostream* ) {}
};
// Output streams // Output streams
static IO::ParallelStreamBuffer DebugStreamBuffer1; #if 0
static IO::ParallelStreamBuffer DebugStreamBuffer2; static IO::ParallelStreamBuffer DebugStreamBuffer1;
static IO::ParallelStreamBuffer DebugStreamBuffer2;
#else
static NullBufferClass DebugStreamBuffer1;
static NullBufferClass DebugStreamBuffer2;
#endif
std::ostream DebugStream1(&DebugStreamBuffer1); std::ostream DebugStream1(&DebugStreamBuffer1);
std::ostream DebugStream2(&DebugStreamBuffer2); std::ostream DebugStream2(&DebugStreamBuffer2);

2
visit/vtkHelpers.h
View File

@ -196,7 +196,7 @@ inline vtkDataSet* meshToVTK( std::shared_ptr<const IO::Mesh> mesh )
mesh2 = DomainToVTK( std::dynamic_pointer_cast<const IO::DomainMesh>(mesh) ); mesh2 = DomainToVTK( std::dynamic_pointer_cast<const IO::DomainMesh>(mesh) );
DebugStream1 << " Volume mesh created" << std::endl; DebugStream1 << " Volume mesh created" << std::endl;
} else { } else {
DebugStream1 << " Error, unknown mesh type" << std::endl; //DebugStream1 << " Error, unknown mesh type" << std::endl;
return NULL; return NULL;
} }
return mesh2; return mesh2;