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

Adding initial code to map bubble ids between timesteps

Browse Source
This commit is contained in:
Mark Berrill 2015-07-16 22:08:43 -04:00
parent 5d86abf774
commit f4a2f22f57
4 changed files with 388 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
IO/Mesh.cpp
View File

@ -326,6 +326,7 @@ std::pair<size_t,void*> DomainMesh::pack( int level ) const
std::pair<size_t,void*> data(0,NULL);
data.first = 7*sizeof(double);
data.second = new double[7];
memset(data.second,0,7*sizeof(double));
int *data_int = reinterpret_cast<int*>(data.second);
double *data_double = &reinterpret_cast<double*>(data.second)[4];
data_int[0] = nprocx;

154
analysis/analysis.cpp
View File

@ -736,4 +736,158 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, RankInfoStruct rank_inf
}
/******************************************************************
* Compute the mapping of blob ids between timesteps *
******************************************************************/
void gatherSet( std::set<int>& set )
{
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);
MPI_Allgather(&send_count,1,MPI_INT,getPtr(recv_count),1,MPI_INT,MPI_COMM_WORLD);
for (int i=1; i<nprocs; i++)
recv_disp[i] = recv_disp[i-1] + recv_count[i-1];
std::vector<int> recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]);
MPI_Allgatherv(getPtr(send_data),send_count,MPI_INT,
getPtr(recv_data),getPtr(recv_count),getPtr(recv_disp),MPI_INT,
MPI_COMM_WORLD);
for (size_t i=0; i<recv_data.size(); i++)
set.insert(recv_data[i]);
}
void gatherSrcIDMap( std::map<int,std::set<int> >& src_map )
{
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;
const std::set<int>& src_ids = it->second;
send_data.push_back(id);
send_data.push_back(src_ids.size());
for (std::set<int>::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2)
send_data.push_back(*it2);
}
int send_count = send_data.size();
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);
for (int i=1; i<nprocs; i++)
recv_disp[i] = recv_disp[i-1] + recv_count[i-1];
std::vector<int> recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]);
MPI_Allgatherv(getPtr(send_data),send_count,MPI_INT,
getPtr(recv_data),getPtr(recv_count),getPtr(recv_disp),MPI_INT,
MPI_COMM_WORLD);
size_t i=0;
while ( i < recv_data.size() ) {
int id = recv_data[i];
int count = recv_data[i+1];
i += 2;
std::set<int>& src_ids = src_map[id];
for (int j=0; j<count; j++,i++)
src_ids.insert(recv_data[i]);
}
}
void addSrcDstIDs( int src_id, std::map<int,std::set<int> >& src_map,
std::map<int,std::set<int> >& dst_map, std::set<int>& src, std::set<int>& dst )
{
src.insert(src_id);
const std::set<int>& dst_ids = dst_map[src_id];
for (std::set<int>::const_iterator it=dst_ids.begin(); it!=dst_ids.end(); ++it) {
if ( dst.find(*it)==dst.end() )
addSrcDstIDs(*it,dst_map,src_map,dst,src);
}
}
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));
PROFILE_START("computeIDMap");
// 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::map<int,std::set<int> > src_map; // Map of the src ids for each dst id
for (size_t i=0; i<ID1.length(); i++) {
if ( ID1(i)>=0 )
src_set.insert(ID1(i));
if ( ID2(i)>=0 )
dst_set.insert(ID2(i));
if ( ID2(i)>=0 && ID1(i)>=0 ) {
std::set<int>& src_ids = src_map[ID2(i)];
src_ids.insert(ID1(i));
}
}
// Communicate the src/dst ids and src id map to all processors and reduce
gatherSet( src_set );
gatherSet( dst_set );
gatherSrcIDMap( src_map );
// Compute the dst id map
std::map<int,std::set<int> > 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) {
int id = it->first;
const std::set<int>& src_ids = it->second;
for (std::set<int>::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2) {
std::set<int>& dst_ids = dst_map[*it2];
dst_ids.insert(id);
}
}
// Perform the mapping of ids
ID_map_struct id_map;
// Find new blobs
for (std::set<int>::const_iterator it=dst_set.begin(); it!=dst_set.end(); ++it) {
if ( src_map.find(*it)==src_map.end() )
id_map.created.push_back(*it);
}
// Fine blobs that disappeared
for (std::set<int>::const_iterator it=src_set.begin(); it!=src_set.end(); ++it) {
if ( dst_map.find(*it)==dst_map.end() )
id_map.destroyed.push_back(*it);
}
// Find blobs with a 1-to-1 mapping
std::vector<int> dst_list;
dst_list.reserve(src_map.size());
for (std::map<int,std::set<int> >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it)
dst_list.push_back(it->first);
for (size_t i=0; i<dst_list.size(); i++) {
int dst_id = dst_list[i];
const std::set<int>& src_ids = src_map[dst_id];
if ( src_ids.size()==1 ) {
int src_id = *src_ids.begin();
const std::set<int>& dst_ids = dst_map[src_id];
if ( dst_ids.size()==1 ) {
ASSERT(*dst_ids.begin()==dst_id);
src_map.erase(dst_id);
dst_map.erase(src_id);
id_map.src_dst.push_back(std::pair<int,int>(src_id,dst_id));
}
}
}
// Handle merge/splits
while ( !dst_map.empty() ) {
// Get a lit of the src-dst ids
std::set<int> 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)
dst_map.erase(*it);
for (std::set<int>::const_iterator it=dst.begin(); it!=dst.end(); ++it)
src_map.erase(*it);
if ( src.size()==1 ) {
// Bubble split
id_map.split.push_back( BlobIDSplitStruct(*src.begin(),std::vector<int>(dst.begin(),dst.end())) );
} else if ( dst.size()==1 ) {
// Bubble merge
id_map.merge.push_back( BlobIDMergeStruct(std::vector<int>(src.begin(),src.end()),*dst.begin()) );
} else {
// Bubble split/merge
id_map.merge_split.push_back( BlobIDMergeSplitStruct(
std::vector<int>(src.begin(),src.end()), std::vector<int>(dst.begin(),dst.end() ) ) );
}
}
ASSERT(src_map.empty());
PROFILE_STOP("computeIDMap");
return id_map;
}

33
analysis/analysis.h
View File

@ -4,6 +4,10 @@
#include "common/Array.h"
#include "common/Communication.h"
#include <set>
#include <map>
#include <vector>
/*!
* @brief Compute the blob
@ -25,6 +29,7 @@ int ComputeBlob( IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator,
const DoubleArray &F, const DoubleArray &S, double vf, double vs, int startx, int starty,
int startz, IntArray &temp, bool periodic=true );
/*!
* @brief Compute the blob
* @details Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
@ -49,10 +54,10 @@ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
* @param[out] ComponentLabel
* @param[in] periodic
*/
int ComputeLocalPhaseComponent(IntArray &PhaseID, int VALUE, IntArray &ComponentLabel,
bool periodic );
/*!
* @brief Compute the blob
* @details Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
@ -72,6 +77,7 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, RankInfoStruct rank_info,
const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS,
IntArray& GlobalBlobID );
/*!
* @brief Compute component of the specified phase
* @details Compute component of specified phase PhaseID=VALUE
@ -85,10 +91,10 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, RankInfoStruct rank_info,
* @param[out] GlobalBlobID The ids of the blobs for the phase
* @return Return the number of components in the specified phase
*/
int ComputeGlobalPhaseComponent( int nx, int ny, int nz, RankInfoStruct rank_info,
IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID );
/*!
* @brief Reorder the blobs
* @details Reorder the blobs based on the number of cells they contain
@ -101,5 +107,28 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, RankInfoStruct rank_inf
void ReorderBlobIDs( IntArray& ID );
typedef std::pair<int,std::vector<int> > BlobIDSplitStruct;
typedef std::pair<std::vector<int>,int> BlobIDMergeStruct;
typedef std::pair<std::vector<int>,std::vector<int> > BlobIDMergeSplitStruct;
struct ID_map_struct {
std::vector<int> created; // list of new blobs that were created
std::vector<int> 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<BlobIDSplitStruct> split; // list of blobs that split
std::vector<BlobIDMergeStruct> merge; // list of blobs that merged
std::vector<BlobIDMergeSplitStruct> merge_split; // list of blobs that both merged and split
};
/*!
* @brief Get the mapping of blob ids between iterations
* @details This functions computes the map of blob ids between iterations
* @return Returns the map of the blob ids. Each final blob may have no source
* ids, one parent, or multiple parents. Each src id may be a parent for multiple blobs.
* @param[in] ID1 The blob ids at the first timestep
* @param[in] ID2 The blob ids at the second timestep
*/
ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 );
#endif

223
tests/TestBlobIdentify.cpp
View File

@ -26,14 +26,14 @@ struct bubble_struct {
Point center;
double radius;
// Get the distance to the bubble
double dist( const Point& p, double Lx, double Ly, double Lz ) {
double dist( const Point& p, double Lx, double Ly, double Lz ) const {
double x = std::min(fabs(p.x-center.x),std::min(fabs(p.x-center.x-Lx),fabs(p.x-center.x+Lx)));
double y = std::min(fabs(p.y-center.y),std::min(fabs(p.y-center.y-Ly),fabs(p.y-center.y+Ly)));
double z = std::min(fabs(p.z-center.z),std::min(fabs(p.z-center.z-Lz),fabs(p.z-center.z+Lz)));
return sqrt(x*x+y*y+z*z)-radius;
}
// Check if this bubble overlaps with rhs
bool overlap( const bubble_struct& rhs, double Lx, double Ly, double Lz ) {
bool overlap( const bubble_struct& rhs, double Lx, double Ly, double Lz ) const {
return dist(rhs.center,Lx,Ly,Lz) <= radius+rhs.radius;
}
// Create a random bubble
@ -73,6 +73,60 @@ std::vector<bubble_struct> create_bubbles( int N_bubbles, double Lx, double Ly,
}
// Fill the Phase/SignDist info from the bubble list
void fillBubbleData( const std::vector<bubble_struct>& bubbles, DoubleArray& Phase,
DoubleArray& SignDist, double Lx, double Ly, double Lz, const RankInfoStruct rank_info )
{
PROFILE_START("fillBubbleData");
int nx = Phase.size(0)-2;
int ny = Phase.size(1)-2;
int nz = Phase.size(2)-2;
Phase.fill(1);
SignDist.fill(0);
for (int k=0; k<nz; k++) {
double z = Lz*(nz*rank_info.kz+k+0.5)/(nz*rank_info.nz);
for (int j=0; j<ny; j++) {
double y = Ly*(ny*rank_info.jy+j+0.5)/(ny*rank_info.ny);
for (int i=0; i<nx; i++) {
double x = Lx*(nx*rank_info.ix+i+0.5)/(nx*rank_info.nx);
double dist = 1e100;
for (size_t b=0; b<bubbles.size(); b++)
dist = std::min(dist,bubbles[b].dist(Point(x,y,z),Lx,Ly,Lz));
SignDist(i+1,j+1,k+1) = -dist;
}
}
}
PROFILE_STOP("fillBubbleData");
}
// Shift all of the data by the given number of cells
void shift_data( DoubleArray& data, int sx, int sy, int sz, const RankInfoStruct& rank_info )
{
int nx = data.size(0)-2;
int ny = data.size(1)-2;
int nz = data.size(2)-2;
int ngx = nx+2*abs(sx);
int ngy = ny+2*abs(sy);
int ngz = nz+2*abs(sz);
Array<double> tmp1(nx,ny,nz), tmp2(ngx,ngy,ngz), tmp3(ngx,ngy,ngz);
fillHalo<double> fillData1(rank_info,nx,ny,nz,1,1,1,0,1);
fillHalo<double> fillData2(rank_info,nx,ny,nz,abs(sx),abs(sy),abs(sz),0,1);
fillData1.copy(data,tmp1);
fillData2.copy(tmp1,tmp2);
fillData2.fill(tmp2);
for (int k=abs(sz); k<nz+abs(sz); k++) {
for (int j=abs(sy); j<ny+abs(sy); j++) {
for (int i=abs(sx); i<nx+abs(sx); i++)
tmp3(i,j,k) = tmp2(i-sx,j-sy,k-sz);
}
}
fillData2.copy(tmp3,tmp1);
fillData1.copy(tmp1,data);
fillData1.fill(data);
}
// Main
int main(int argc, char **argv)
{
@ -81,12 +135,10 @@ int main(int argc, char **argv)
MPI_Init(&argc,&argv);
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
#ifdef USE_TIMER
PROFILE_ENABLE(0);
PROFILE_DISABLE_TRACE();
PROFILE_SYNCHRONIZE();
PROFILE_START("main");
#endif
if ( rank==0 ) {
printf("-----------------------------------------------------------\n");
printf("Testing Blob Identification \n");
@ -109,22 +161,8 @@ int main(int argc, char **argv)
// Create the dummy info
DoubleArray Phase(nx+2,ny+2,nz+2);
DoubleArray SignDist(nx+2,ny+2,nz+2);
Phase.fill(1);
SignDist(0);
std::vector<bubble_struct> bubbles = create_bubbles(20,Lx,Ly,Lz);
for (int k=0; k<nz; k++) {
double z = Lz*(nz*rank_info.kz+k+0.5)/(nz*nproc[2]);
for (int j=0; j<ny; j++) {
double y = Ly*(ny*rank_info.jy+j+0.5)/(ny*nproc[1]);
for (int i=0; i<nx; i++) {
double x = Lx*(nx*rank_info.ix+i+0.5)/(nx*nproc[0]);
double dist = 1e100;
for (size_t b=0; b<bubbles.size(); b++)
dist = std::min(dist,bubbles[b].dist(Point(x,y,z),Lx,Ly,Lz));
SignDist(i+1,j+1,k+1) = -dist;
}
}
}
fillBubbleData( bubbles, Phase, SignDist, Lx, Ly, Lz, rank_info );
// Communication the halos
fillHalo<double> fillData(rank_info,nx,ny,nz,1,1,1,0,1);
@ -168,6 +206,7 @@ int main(int argc, char **argv)
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( 0, meshData, 2 );
int save_it = 1;
// Check the results
int N_errors = 0;
@ -176,13 +215,153 @@ int main(int argc, char **argv)
N_errors++;
}
// Move the blobs and connect them to the previous ids
PROFILE_START("constant velocity test");
if ( rank==0 ) { printf("Running constant velocity blob test\n"); }
for (int i=0; i<20; i++, save_it++) {
// Shift all the data
shift_data( Phase, 3, -2, 1, rank_info );
shift_data( SignDist, 3, -2, 1, rank_info );
// Find blob domains
IntArray GlobalBlobID2;
int nblobs2 = ComputeGlobalBlobIDs(nx,ny,nz,rank_info,
Phase,SignDist,vF,vS,GlobalBlobID2);
if ( nblobs2 != nblobs ) {
printf("Error, number of blobs changed under constant velocity (%i,%i)\n",nblobs,nblobs2);
N_errors++;
}
// Identify the blob maps and renumber the ids
ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2);
bool pass = (int)map.src_dst.size()==nblobs;
pass = pass && map.created.empty();
pass = pass && map.destroyed.empty();
for (size_t j=0; j<map.src_dst.size(); j++)
pass = pass && map.src_dst[j].first==map.src_dst[j].second;
pass = pass && map.split.empty();
pass = pass && map.merge.empty();
pass = pass && map.merge_split.empty();
if ( !pass ) {
printf("Error, blob ids do not match in constant velocity test\n");
N_errors++;
}
GlobalBlobID = GlobalBlobID2;
// Save the results
fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( save_it, meshData, 2 );
}
PROFILE_STOP("constant velocity test");
// Move the bubbles in independent directions to test merge/splitting of bubbles
PROFILE_START("moving bubble test");
if ( rank==0 ) { printf("Running moving bubble test\n"); }
std::vector<Point> velocity(bubbles.size());
if ( rank==0 ) {
for (size_t i=0; i<bubbles.size(); i++) {
velocity[i].x = bubbles[i].radius*(2*rand2()-1);
velocity[i].y = bubbles[i].radius*(2*rand2()-1);
velocity[i].z = bubbles[i].radius*(2*rand2()-1);
}
}
MPI_Bcast(&velocity[0],bubbles.size()*sizeof(Point),MPI_CHAR,0,MPI_COMM_WORLD);
fillBubbleData( bubbles, Phase, SignDist, Lx, Ly, Lz, rank_info );
fillData.fill(Phase);
fillData.fill(SignDist);
ComputeGlobalBlobIDs(nx,ny,nz,rank_info,Phase,SignDist,vF,vS,GlobalBlobID);
fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( save_it, meshData, 2 );
save_it++;
for (int i=0; i<25; i++, save_it++) {
// Move the bubbles
for (size_t j=0; j<bubbles.size(); j++) {
bubbles[j].center.x = fmod(bubbles[j].center.x+velocity[j].x+Lx,Lx);
bubbles[j].center.y = fmod(bubbles[j].center.y+velocity[j].y+Ly,Ly);
bubbles[j].center.z = fmod(bubbles[j].center.z+velocity[j].z+Lz,Lz);
}
fillBubbleData( bubbles, Phase, SignDist, Lx, Ly, Lz, rank_info );
fillData.fill(Phase);
fillData.fill(SignDist);
// Compute the ids
IntArray GlobalBlobID2;
int nblobs2 = ComputeGlobalBlobIDs(nx,ny,nz,rank_info,Phase,SignDist,vF,vS,GlobalBlobID2);
// Identify the blob maps and renumber the ids
ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2);
int N1 = 0;
int N2 = 0;
if ( rank==0 ) {
if ( !map.destroyed.empty() ) {
N1 += map.destroyed.size();
printf(" %i: destroyed:",i+1);
for (size_t j=0; j<map.destroyed.size(); j++)
printf(" %i",map.destroyed[j]);
printf("\n");
}
if ( !map.created.empty() ) {
N2 += map.created.size();
printf(" %i: created:",i+1);
for (size_t j=0; j<map.created.size(); j++)
printf(" %i",map.created[j]);
printf("\n");
}
N1 += map.src_dst.size();
N2 += map.src_dst.size();
for (size_t j=0; j<map.split.size(); j++ ) {
N1 += 1;
N2 += map.split[j].second.size();
printf(" %i: split: %i -",i+1,map.split[j].first);
for (size_t k=0; k<map.split[j].second.size(); k++)
printf(" %i",map.split[j].second[k]);
printf("\n");
}
for (size_t j=0; j<map.merge.size(); j++ ) {
N1 += map.merge[j].first.size();
N2 += 1;
printf(" %i: merged:",i+1);
for (size_t k=0; k<map.merge[j].first.size(); k++)
printf(" %i",map.merge[j].first[k]);
printf(" - %i\n",map.merge[j].second);
}
for (size_t j=0; j<map.merge_split.size(); j++ ) {
N1 += map.merge_split[j].first.size();
N2 += map.merge_split[j].second.size();
printf(" %i: merged/split:",i+1);
for (size_t k=0; k<map.merge_split[j].first.size(); k++)
printf(" %i",map.merge_split[j].first[k]);
printf(" -");
for (size_t k=0; k<map.merge_split[j].second.size(); k++)
printf(" %i",map.merge_split[j].second[k]);
printf("\n");
}
}
MPI_Bcast(&N1,1,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(&N2,1,MPI_INT,0,MPI_COMM_WORLD);
if ( N1!=nblobs || N2!=nblobs2 ) {
if ( rank==0 )
printf("Error, blob ids do not map in moving bubble test (%i,%i,%i,%i)\n",
nblobs,nblobs2,N1,N2);
N_errors++;
}
nblobs = nblobs2;
GlobalBlobID = GlobalBlobID2;
// Save the results
fillData.copy(Phase,PhaseVar->data);
fillData.copy(SignDist,SignDistVar->data);
fillData.copy(GlobalBlobID,BlobIDVar->data);
IO::writeData( save_it, meshData, 2 );
}
PROFILE_STOP("moving bubble test");
// Finished
#ifdef USE_TIMER
PROFILE_STOP("main");
PROFILE_SAVE("TestBlobIdentify",false);
#endif
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
return 0;
return N_errors;
}