Adding parallel blob identification
This commit is contained in:
parent
ef25839769
commit
b7634dc099
|
@ -124,6 +124,7 @@ ENDMACRO()
|
|||
IF ( NOT ONLY_BUILD_DOCS )
|
||||
BEGIN_PACKAGE_CONFIG( lbpm-wia )
|
||||
ADD_PACKAGE_SUBDIRECTORY( common )
|
||||
ADD_PACKAGE_SUBDIRECTORY( analysis )
|
||||
ADD_PACKAGE_SUBDIRECTORY( IO )
|
||||
IF ( USE_CUDA )
|
||||
ADD_PACKAGE_SUBDIRECTORY( gpu )
|
||||
|
|
488
analysis/analysis.cpp
Normal file
488
analysis/analysis.cpp
Normal file
|
@ -0,0 +1,488 @@
|
|||
#include "analysis/analysis.h"
|
||||
#include "ProfilerApp.h"
|
||||
#include <iostream>
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
inline TYPE* getPtr( std::vector<TYPE>& x ) { return x.empty() ? NULL:&x[0]; }
|
||||
template<class TYPE>
|
||||
inline const TYPE* getPtr( const std::vector<TYPE>& x ) { return x.empty() ? NULL:&x[0]; }
|
||||
|
||||
|
||||
/******************************************************************
|
||||
* Compute the blobls *
|
||||
******************************************************************/
|
||||
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)
|
||||
{
|
||||
// Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
|
||||
// F>vf => oil phase S>vs => in porespace
|
||||
// update the list of blobs, indicator mesh
|
||||
int m = F.size(0); // maxima for the meshes
|
||||
int n = F.size(1);
|
||||
int o = F.size(2);
|
||||
|
||||
int cubes_in_blob=0;
|
||||
int nrecent = 1; // number of nodes added at most recent sweep
|
||||
temp(0,0) = startx; // Set the initial point as a "seed" for the sweeps
|
||||
temp(1,0) = starty;
|
||||
temp(2,0) = startz;
|
||||
int ntotal = 1; // total number of nodes in blob
|
||||
indicator(startx,starty,startz) = nblobs;
|
||||
|
||||
int p,s,x,y,z,start,finish,nodx,nody,nodz;
|
||||
int imin=startx,imax=startx,jmin=starty,jmax=starty; // initialize maxima / minima
|
||||
int kmin=startz,kmax=startz;
|
||||
int d[26][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},
|
||||
{1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0},{1,0,1},{-1,0,1},
|
||||
{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-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}}; // directions to neighbors
|
||||
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
|
||||
bool status = 1; // status == true => continue to look for points
|
||||
while (status == 1){
|
||||
start = ntotal - nrecent;
|
||||
finish = ntotal;
|
||||
nrecent = 0; // set recent points back to zero for next sweep through
|
||||
for (s=start;s<finish;s++){
|
||||
// Loop over recent points; look for new points
|
||||
x = temp(0,s);
|
||||
y = temp(1,s);
|
||||
z = temp(2,s);
|
||||
// Looop over the directions
|
||||
for (p=0;p<26;p++){
|
||||
nodx=x+d[p][0];
|
||||
nody=y+d[p][1];
|
||||
nodz=z+d[p][2];
|
||||
if ( periodic ) {
|
||||
if (nodx < 0 ){ nodx = m-1; } // Periodic BC for x
|
||||
if (nodx > m-1 ){ nodx = 0; }
|
||||
if (nody < 0 ){ nody = n-1; } // Periodic BC for y
|
||||
if (nody > n-1 ){ nody = 0; }
|
||||
if (nodz < 0 ){ nodz = o-1; } // Periodic BC for z
|
||||
if (nodz > o-1 ){ nodz = 0; }
|
||||
} else {
|
||||
if ( nodx<0 || nodx>=m || nody<0 || nody>=n || nodz<0 || nodz>=o )
|
||||
continue;
|
||||
}
|
||||
if ( F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs
|
||||
&& indicator(nodx,nody,nodz) == -1 ){
|
||||
// Node is a part of the blob - add it to the list
|
||||
temp(0,ntotal) = nodx;
|
||||
temp(1,ntotal) = nody;
|
||||
temp(2,ntotal) = nodz;
|
||||
ntotal++;
|
||||
nrecent++;
|
||||
// Update the indicator map
|
||||
indicator(nodx,nody,nodz) = nblobs;
|
||||
// Update the min / max for the cube loop
|
||||
if ( nodx < imin ){ imin = nodx; }
|
||||
if ( nodx > imax ){ imax = nodx; }
|
||||
if ( nody < jmin ){ jmin = nody; }
|
||||
if ( nody > jmax ){ jmax = nody; }
|
||||
if ( nodz < kmin ){ kmin = nodz; }
|
||||
if ( nodz > kmax ){ kmax = nodz; }
|
||||
}
|
||||
else if (F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs
|
||||
&& indicator(nodx,nody,nodz) > -1 && indicator(nodx,nody,nodz) != nblobs){
|
||||
// Some kind of error in algorithm
|
||||
printf("Error in blob search algorithm!");
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
if ( nrecent == 0){
|
||||
status = 0;
|
||||
}
|
||||
}
|
||||
// Use points in temporary storage array to add cubes to the list of blobs
|
||||
if ( imin > 0) { imin = imin-1; }
|
||||
// if ( imax < m-1) { imax = imax+1; }
|
||||
if ( jmin > 0) { jmin = jmin-1; }
|
||||
// if ( jmax < n-1) { jmax = jmax+1; }
|
||||
if ( kmin > 0) { kmin = kmin-1; }
|
||||
// if ( kmax < o-1) { kmax = kmax+1; }
|
||||
int i,j,k;
|
||||
bool add;
|
||||
for (k=kmin;k<kmax;k++){
|
||||
for (j=jmin;j<jmax;j++){
|
||||
for (i=imin;i<imax;i++){
|
||||
// If cube(i,j,k) has any nodes in blob, add it to the list
|
||||
// Loop over cube edges
|
||||
add = 0;
|
||||
for (p=0;p<8;p++){
|
||||
nodx = i+cube[p][0];
|
||||
nody = j+cube[p][1];
|
||||
nodz = k+cube[p][2];
|
||||
if ( indicator(nodx,nody,nodz) == nblobs ){
|
||||
// Cube corner is in this blob
|
||||
add = 1;
|
||||
}
|
||||
}
|
||||
if (add == 1){
|
||||
// Add cube to the list
|
||||
blobs(0,ncubes) = i;
|
||||
blobs(1,ncubes) = j;
|
||||
blobs(2,ncubes) = k;
|
||||
ncubes++;
|
||||
cubes_in_blob++;
|
||||
// Loop again to check for overlap
|
||||
for (p=0;p<8;p++){
|
||||
nodx = i+cube[p][0];
|
||||
nody = j+cube[p][1];
|
||||
nodz = k+cube[p][2];
|
||||
if (indicator(nodx,nody,nodz) > -1 && indicator(nodx,nody,nodz) != nblobs){
|
||||
printf("Overlapping cube!");
|
||||
std::cout << i << ", " << j << ", " << k << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return cubes_in_blob;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/******************************************************************
|
||||
* Compute the local blob ids *
|
||||
******************************************************************/
|
||||
int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
|
||||
double vF, double vS, IntArray& LocalBlobID, bool periodic )
|
||||
{
|
||||
PROFILE_START("ComputeLocalBlobIDs");
|
||||
size_t Nx = Phase.size(0);
|
||||
size_t Ny = Phase.size(1);
|
||||
size_t Nz = Phase.size(2);
|
||||
ASSERT(SignDist.size(0)==Nx&&SignDist.size(1)==Ny&&SignDist.size(2)==Nz);
|
||||
// Initialize output
|
||||
LocalBlobID.resize(Nx,Ny,Nz);
|
||||
// Compute the local blob ids
|
||||
const int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
|
||||
size_t N = Nx*Ny*Nz;
|
||||
int nblobs = 0;
|
||||
int ncubes = 0; // total number of nodes in any blob
|
||||
IntArray blobs(3,N); // store indices for blobs (cubes)
|
||||
IntArray temp(3,N); // temporary storage array
|
||||
IntArray b(N); // number of nodes in each blob
|
||||
for (size_t k=0; k<Nz; k++ ) {
|
||||
for (size_t j=0; j<Ny; j++) {
|
||||
for (size_t i=0; i<Nx; i++) {
|
||||
if ( SignDist(i,j,k) < 0.0) {
|
||||
// Solid phase
|
||||
LocalBlobID(i,j,k) = -2;
|
||||
} else{
|
||||
LocalBlobID(i,j,k) = -1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
for (size_t k=0; k<Nz; k++ ) {
|
||||
for (size_t j=0; j<Ny; j++) {
|
||||
for (size_t i=0; i<Nx; i++) {
|
||||
if ( LocalBlobID(i,j,k)==-1 && Phase(i,j,k)>vF && SignDist(i,j,k)>vS ) {
|
||||
// node i,j,k is in the porespace
|
||||
b(nblobs) = ComputeBlob(blobs,nblobs,ncubes,LocalBlobID,
|
||||
Phase,SignDist,vF,vS,i,j,k,temp,periodic);
|
||||
nblobs++;
|
||||
}
|
||||
if ( nblobs > (int)b.length()-1){
|
||||
printf("Increasing size of blob list \n");
|
||||
b.resize(2*b.length());
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Go over all cubes again -> add any that do not contain nw phase
|
||||
bool add=1; // Set to false if any corners contain nw-phase ( F > vF)
|
||||
size_t count_in=0,count_out=0;
|
||||
size_t nodx,nody,nodz;
|
||||
for (size_t k=0; k<Nz-1; k++ ) {
|
||||
for (size_t j=0; j<Ny-1; j++) {
|
||||
for (size_t i=0; i<Nx-1; i++) {
|
||||
// Loop over cube corners
|
||||
int add=1; // initialize to true - add unless corner occupied by nw-phase
|
||||
for (int p=0; p<8; p++) {
|
||||
nodx=i+cube[p][0];
|
||||
nody=j+cube[p][1];
|
||||
nodz=k+cube[p][2];
|
||||
if ( LocalBlobID(nodx,nody,nodz) > -1 ){
|
||||
// corner occupied by nw-phase -> do not add
|
||||
add = 0;
|
||||
}
|
||||
}
|
||||
if ( add == 1 ){
|
||||
blobs(0,ncubes) = i;
|
||||
blobs(1,ncubes) = j;
|
||||
blobs(2,ncubes) = k;
|
||||
ncubes++;
|
||||
count_in++;
|
||||
}
|
||||
else { count_out++; }
|
||||
}
|
||||
}
|
||||
}
|
||||
b(nblobs) = count_in;
|
||||
PROFILE_STOP("ComputeLocalBlobIDs");
|
||||
return nblobs;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/******************************************************************
|
||||
* Reorder the global blob ids *
|
||||
******************************************************************/
|
||||
static int ReorderBlobIDs2( IntArray& ID, int N_blobs, int ngx, int ngy, int ngz )
|
||||
{
|
||||
if ( N_blobs==0 )
|
||||
return 0;
|
||||
PROFILE_START("ReorderBlobIDs2",1);
|
||||
ASSERT(sizeof(long long int)==sizeof(int64_t));
|
||||
double *local_size = new double[N_blobs];
|
||||
double *global_size = new double[N_blobs];
|
||||
for (int i=0; i<N_blobs; i++)
|
||||
local_size[i] = 0;
|
||||
for (int i=0; i<N_blobs; i++)
|
||||
global_size[i] = 0;
|
||||
int max_id = -1;
|
||||
for (size_t k=ngz; k<ID.size(2)-ngz; k++) {
|
||||
for (size_t j=ngy; j<ID.size(1)-ngy; j++) {
|
||||
for (size_t i=ngx; i<ID.size(0)-ngx; i++) {
|
||||
int id = ID(i,j,k);
|
||||
if ( id >= 0 )
|
||||
local_size[id] += 1;
|
||||
max_id = std::max(max_id,id);
|
||||
}
|
||||
}
|
||||
}
|
||||
ASSERT(max_id<N_blobs);
|
||||
MPI_Allreduce(local_size,global_size,N_blobs,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
|
||||
std::vector<std::pair<double,int> > map1(N_blobs);
|
||||
int N_blobs2 = 0;
|
||||
for (int i=0; i<N_blobs; i++) {
|
||||
map1[i].first = global_size[i];
|
||||
map1[i].second = i;
|
||||
if ( global_size[i] > 0 )
|
||||
N_blobs2++;
|
||||
}
|
||||
std::sort( map1.begin(), map1.end() );
|
||||
std::vector<int> map2(N_blobs,-1);
|
||||
for (int i=0; i<N_blobs; i++) {
|
||||
map2[map1[N_blobs-i-1].second] = i;
|
||||
}
|
||||
for (size_t i=0; i<ID.length(); i++) {
|
||||
if ( ID(i) >= 0 )
|
||||
ID(i) = map2[ID(i)];
|
||||
}
|
||||
delete [] local_size;
|
||||
delete [] global_size;
|
||||
PROFILE_STOP("ReorderBlobIDs2",1);
|
||||
return N_blobs2;
|
||||
}
|
||||
void ReorderBlobIDs( IntArray& ID )
|
||||
{
|
||||
PROFILE_START("ReorderBlobIDs");
|
||||
int tmp = ID.max()+1;
|
||||
int N_blobs = 0;
|
||||
MPI_Allreduce(&tmp,&N_blobs,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
|
||||
ReorderBlobIDs2(ID,N_blobs,1,1,1);
|
||||
PROFILE_STOP("ReorderBlobIDs");
|
||||
}
|
||||
|
||||
|
||||
|
||||
/******************************************************************
|
||||
* Compute the global blob ids *
|
||||
******************************************************************/
|
||||
struct global_id_info_struct {
|
||||
int64_t new_id;
|
||||
std::set<int64_t> remote_ids;
|
||||
};
|
||||
// Send the local ids and their new value to all neighbors
|
||||
static void updateRemoteIds(
|
||||
const std::map<int64_t,global_id_info_struct>& map,
|
||||
const std::vector<int>& neighbors,
|
||||
int N_send, const std::vector<int>& N_recv,
|
||||
int64_t *send_buf, std::vector<int64_t*>& recv_buf,
|
||||
std::map<int64_t,int64_t>& remote_map )
|
||||
{
|
||||
std::vector<MPI_Request> send_req(neighbors.size());
|
||||
std::vector<MPI_Request> recv_req(neighbors.size());
|
||||
std::vector<MPI_Status> status(neighbors.size());
|
||||
std::map<int64_t,global_id_info_struct>::const_iterator it = map.begin();
|
||||
ASSERT(N_send==(int)map.size());
|
||||
for (size_t i=0; i<map.size(); i++, ++it) {
|
||||
send_buf[2*i+0] = it->first;
|
||||
send_buf[2*i+1] = it->second.new_id;
|
||||
}
|
||||
for (size_t i=0; i<neighbors.size(); i++) {
|
||||
MPI_Isend( send_buf, 2*N_send, MPI_LONG_LONG, neighbors[i], 0, MPI_COMM_WORLD, &send_req[i] );
|
||||
MPI_Irecv( recv_buf[i], 2*N_recv[i], MPI_LONG_LONG, neighbors[i], 0, MPI_COMM_WORLD, &recv_req[i] );
|
||||
}
|
||||
for (it=map.begin(); it!=map.end(); ++it) {
|
||||
remote_map[it->first] = it->second.new_id;
|
||||
}
|
||||
for (size_t i=0; i<neighbors.size(); i++) {
|
||||
MPI_Wait(&recv_req[i],&status[i]);
|
||||
for (int j=0; j<N_recv[i]; j++)
|
||||
remote_map[recv_buf[i][2*j+0]] = recv_buf[i][2*j+1];
|
||||
}
|
||||
MPI_Waitall(neighbors.size(),getPtr(send_req),getPtr(status));
|
||||
}
|
||||
// Compute a new local id for each local id
|
||||
static bool updateLocalIds( const std::map<int64_t,int64_t>& remote_map,
|
||||
std::map<int64_t,global_id_info_struct>& map )
|
||||
{
|
||||
bool changed = false;
|
||||
std::map<int64_t,global_id_info_struct>::iterator it;
|
||||
for (it=map.begin(); it!=map.end(); ++it) {
|
||||
int64_t id = it->second.new_id;
|
||||
std::set<int64_t>::const_iterator it2;
|
||||
for (it2=it->second.remote_ids.begin(); it2!=it->second.remote_ids.end(); ++it2) {
|
||||
int64_t id2 = remote_map.find(*it2)->second;
|
||||
id = std::min(id,id2);
|
||||
}
|
||||
changed = changed || it->second.new_id!=id;
|
||||
it->second.new_id = id;
|
||||
}
|
||||
return changed;
|
||||
}
|
||||
int ComputeGlobalBlobIDs( int nx, int ny, int nz, RankInfoStruct rank_info,
|
||||
const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS,
|
||||
IntArray& GlobalBlobID )
|
||||
{
|
||||
PROFILE_START("ComputeGlobalBlobIDs");
|
||||
const int rank = rank_info.rank[1][1][1];
|
||||
int nprocs;
|
||||
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
|
||||
// First compute the local ids
|
||||
IntArray LocalIDs;
|
||||
int nblobs = ComputeLocalBlobIDs(Phase,SignDist,vF,vS,LocalIDs,false);
|
||||
std::vector<int> N_blobs(nprocs,0);
|
||||
PROFILE_START("ComputeGlobalBlobIDs-Allgather",1);
|
||||
MPI_Allgather(&nblobs,1,MPI_INT,getPtr(N_blobs),1,MPI_INT,MPI_COMM_WORLD);
|
||||
PROFILE_STOP("ComputeGlobalBlobIDs-Allgather",1);
|
||||
int64_t N_blobs_tot = 0;
|
||||
int offset = 0;
|
||||
for (int i=0; i<rank; i++)
|
||||
offset += N_blobs[i];
|
||||
for (int i=0; i<nprocs; i++)
|
||||
N_blobs_tot += N_blobs[i];
|
||||
INSIST(N_blobs_tot<0x80000000,"Maximum number of blobs exceeded");
|
||||
// Compute temporary global ids
|
||||
for (size_t i=0; i<LocalIDs.length(); i++) {
|
||||
if ( LocalIDs(i) >= 0 )
|
||||
LocalIDs(i) += offset;
|
||||
}
|
||||
// Copy the ids and get the neighbors through the halos
|
||||
GlobalBlobID = LocalIDs;
|
||||
fillHalo<int> fillData(rank_info,nx,ny,nz,1,1,1,0,1,true,true,true);
|
||||
fillData.fill(GlobalBlobID);
|
||||
// Create a list of all neighbor ranks (excluding self)
|
||||
std::vector<int> neighbors;
|
||||
neighbors.push_back( rank_info.rank[0][1][1] );
|
||||
neighbors.push_back( rank_info.rank[2][1][1] );
|
||||
neighbors.push_back( rank_info.rank[1][0][1] );
|
||||
neighbors.push_back( rank_info.rank[1][2][1] );
|
||||
neighbors.push_back( rank_info.rank[1][1][0] );
|
||||
neighbors.push_back( rank_info.rank[1][1][2] );
|
||||
std::unique(neighbors.begin(),neighbors.end());
|
||||
if ( std::find(neighbors.begin(),neighbors.end(),rank) != neighbors.end() )
|
||||
neighbors.erase(std::find(neighbors.begin(),neighbors.end(),rank));
|
||||
// Create a map of all local ids to the neighbor ids
|
||||
std::map<int64_t,global_id_info_struct> map;
|
||||
std::set<int64_t> local;
|
||||
for (size_t i=0; i<LocalIDs.length(); i++) {
|
||||
if ( LocalIDs(i)>=0 ) {
|
||||
local.insert(LocalIDs(i));
|
||||
if ( LocalIDs(i)!=GlobalBlobID(i) )
|
||||
map[LocalIDs(i)].remote_ids.insert(GlobalBlobID(i));
|
||||
}
|
||||
}
|
||||
std::map<int64_t,global_id_info_struct>::iterator it;
|
||||
for (it=map.begin(); it!=map.end(); ++it) {
|
||||
it->second.new_id = it->first;
|
||||
local.erase(it->first);
|
||||
}
|
||||
// Get the number of ids we will recieve from each rank
|
||||
int N_send = map.size();
|
||||
std::vector<int> N_recv(neighbors.size(),0);
|
||||
std::vector<MPI_Request> send_req(neighbors.size());
|
||||
std::vector<MPI_Request> recv_req(neighbors.size());
|
||||
std::vector<MPI_Status> status(neighbors.size());
|
||||
for (size_t i=0; i<neighbors.size(); i++) {
|
||||
MPI_Isend( &N_send, 1, MPI_INT, neighbors[i], 0, MPI_COMM_WORLD, &send_req[i] );
|
||||
MPI_Irecv( &N_recv[i], 1, MPI_INT, neighbors[i], 0, MPI_COMM_WORLD, &recv_req[i] );
|
||||
}
|
||||
MPI_Waitall(neighbors.size(),getPtr(send_req),getPtr(status));
|
||||
MPI_Waitall(neighbors.size(),getPtr(recv_req),getPtr(status));
|
||||
// Allocate memory for communication
|
||||
int64_t *send_buf = new int64_t[2*N_send];
|
||||
std::vector<int64_t*> recv_buf(neighbors.size());
|
||||
for (size_t i=0; i<neighbors.size(); i++)
|
||||
recv_buf[i] = new int64_t[2*N_recv[i]];
|
||||
// Compute a map for the remote ids, and new local id for each id
|
||||
std::map<int64_t,int64_t> remote_map;
|
||||
for (it=map.begin(); it!=map.end(); ++it) {
|
||||
int64_t id = it->first;
|
||||
std::set<int64_t>::const_iterator it2;
|
||||
for (it2=it->second.remote_ids.begin(); it2!=it->second.remote_ids.end(); ++it2) {
|
||||
int64_t id2 = *it2;
|
||||
id = std::min(id,id2);
|
||||
remote_map.insert(std::pair<int64_t,int64_t>(id2,id2));
|
||||
}
|
||||
it->second.new_id = id;
|
||||
}
|
||||
// Iterate until we are done
|
||||
int iteration = 1;
|
||||
PROFILE_START("ComputeGlobalBlobIDs-loop",1);
|
||||
while ( 1 ) {
|
||||
iteration++;
|
||||
// Send the local ids and their new value to all neighbors
|
||||
updateRemoteIds( map, neighbors, N_send, N_recv,send_buf, recv_buf, remote_map );
|
||||
// Compute a new local id for each local id
|
||||
bool changed = updateLocalIds( remote_map, map );
|
||||
// Check if we are finished
|
||||
int test = changed ? 1:0;
|
||||
int result = 0;
|
||||
MPI_Allreduce(&test,&result,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
|
||||
if ( result==0 )
|
||||
break;
|
||||
}
|
||||
PROFILE_STOP("ComputeGlobalBlobIDs-loop",1);
|
||||
// Relabel the ids
|
||||
std::vector<int> final_map(nblobs,-1);
|
||||
for (it=map.begin(); it!=map.end(); ++it)
|
||||
final_map[it->first-offset] = it->second.new_id;
|
||||
for (std::set<int64_t>::const_iterator it2=local.begin(); it2!=local.end(); ++it2)
|
||||
final_map[*it2-offset] = *it2;
|
||||
int ngx = (GlobalBlobID.size(0)-nx)/2;
|
||||
int ngy = (GlobalBlobID.size(1)-ny)/2;
|
||||
int ngz = (GlobalBlobID.size(2)-nz)/2;
|
||||
for (size_t k=ngz; k<GlobalBlobID.size(2)-ngz; k++) {
|
||||
for (size_t j=ngy; j<GlobalBlobID.size(1)-ngy; j++) {
|
||||
for (size_t i=ngx; i<GlobalBlobID.size(0)-ngx; i++) {
|
||||
int id = GlobalBlobID(i,j,k);
|
||||
if ( id >= 0 )
|
||||
GlobalBlobID(i,j,k) = final_map[id-offset];
|
||||
}
|
||||
}
|
||||
}
|
||||
// Fill the ghosts
|
||||
fillHalo<int> fillData2(rank_info,nx,ny,nz,1,1,1,0,1,true,true,true);
|
||||
fillData2.fill(GlobalBlobID);
|
||||
// Reorder based on size (and compress the id space
|
||||
int N_blobs_global = ReorderBlobIDs2(GlobalBlobID,N_blobs_tot,ngx,ngy,ngz);
|
||||
// Finished
|
||||
delete [] send_buf;
|
||||
for (size_t i=0; i<neighbors.size(); i++)
|
||||
delete [] recv_buf[i];
|
||||
PROFILE_STOP("ComputeGlobalBlobIDs");
|
||||
return N_blobs_global;
|
||||
}
|
||||
|
||||
|
||||
|
78
analysis/analysis.h
Normal file
78
analysis/analysis.h
Normal file
|
@ -0,0 +1,78 @@
|
|||
#ifndef COMMON_H_INC
|
||||
#define COMMON_H_INC
|
||||
|
||||
#include "common/Array.h"
|
||||
#include "common/Communication.h"
|
||||
|
||||
|
||||
/*!
|
||||
* @brief Compute the blob
|
||||
* @details Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
|
||||
* @return Returns the number of cubes in the blob
|
||||
* @param[out] blobs blobs
|
||||
* @param[out] nblobs Number of blobs
|
||||
* @param[out] ncubes Number of cubes
|
||||
* @param[out] indicator indicator
|
||||
* @param[in] F F
|
||||
* @param[in] S S
|
||||
* @param[in] vf vf
|
||||
* @param[in] vs vs
|
||||
* @param[in] startx startx
|
||||
* @param[in] starty starty
|
||||
* @param[in/out] temp temp
|
||||
*/
|
||||
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
|
||||
* @return Returns the number of cubes in the blob
|
||||
* @param[in] Phase Phase
|
||||
* @param[in] SignDist SignDist
|
||||
* @param[in] vF vF
|
||||
* @param[in] vS vS
|
||||
* @param[in] S S
|
||||
* @param[out] LocalBlobID The ids of the blobs
|
||||
* @return Returns the number of blobs
|
||||
*/
|
||||
int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
|
||||
double vF, double vS, IntArray& LocalBlobID, bool periodic=true );
|
||||
|
||||
|
||||
/*!
|
||||
* @brief Compute the blob
|
||||
* @details Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
|
||||
* @return Returns the number of cubes in the blob
|
||||
* @param[in] nx Number of elements in the x-direction
|
||||
* @param[in] ny Number of elements in the y-direction
|
||||
* @param[in] nz Number of elements in the z-direction
|
||||
* @param[in] Phase Phase
|
||||
* @param[in] SignDist SignDist
|
||||
* @param[in] vF vF
|
||||
* @param[in] vS vS
|
||||
* @param[in] S S
|
||||
* @param[out] LocalBlobID The ids of the blobs
|
||||
* @return Returns the number of blobs
|
||||
*/
|
||||
int ComputeGlobalBlobIDs( int nx, int ny, int nz, RankInfoStruct rank_info,
|
||||
const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS,
|
||||
IntArray& GlobalBlobID );
|
||||
|
||||
|
||||
/*!
|
||||
* @brief Reorder the blobs
|
||||
* @details Reorder the blobs based on the number of cells they contain
|
||||
* largest first.
|
||||
* @param[in] nx Number of elements in the x-direction
|
||||
* @param[in] ny Number of elements in the y-direction
|
||||
* @param[in] nz Number of elements in the z-direction
|
||||
* @param[in/out] ID The ids of the blobs
|
||||
*/
|
||||
void ReorderBlobIDs( IntArray& ID );
|
||||
|
||||
|
||||
|
||||
#endif
|
|
@ -58,10 +58,6 @@ FUNCTION( CONFIGURE_SHARED_PTR INSTALL_DIR NAMESPACE )
|
|||
IF ( NOT NAMESPACE )
|
||||
SET( NAMESPACE " " )
|
||||
ENDIF()
|
||||
SET( BOOST_SHARED_PTR 0 )
|
||||
SET( MEMORY_SHARED_PTR 0 )
|
||||
SET( MEMORY_TR1_SHARED_PTR 0 )
|
||||
SET( TR1_MEMORY_TR1_SHARED_PTR 0 )
|
||||
IF ( BOOST_SHARED_PTR )
|
||||
FILE(WRITE "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include \"boost/shared_ptr.hpp\"\n")
|
||||
FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include \"boost/weak_ptr.hpp\"\n")
|
||||
|
@ -122,19 +118,20 @@ FUNCTION( WRITE_DUMMY_SHARED_PTR NAMESPACE FILENAME )
|
|||
FILE(APPEND "${FILENAME}" "template<class T> void DefaultDeleter(T* p) {delete p;}\n\n")
|
||||
FILE(APPEND "${FILENAME}" "template<class T> class shared_ptr {\n")
|
||||
FILE(APPEND "${FILENAME}" "public:\n")
|
||||
FILE(APPEND "${FILENAME}" " shared_ptr( ): obj(NULL), count(NULL) {}\n")
|
||||
FILE(APPEND "${FILENAME}" " typedef void (*D)(T*);\n")
|
||||
FILE(APPEND "${FILENAME}" " shared_ptr( ): obj(NULL), deleter(DefaultDeleter<T>), count(NULL) {}\n")
|
||||
FILE(APPEND "${FILENAME}" " shared_ptr( T *ptr, void (*D)(T*)=DefaultDeleter<T>):\n")
|
||||
FILE(APPEND "${FILENAME}" " obj(ptr), deleter(D), count(NULL) { if (ptr) { count = new int; (*count)=1; } } \n")
|
||||
FILE(APPEND "${FILENAME}" " shared_ptr( const shared_ptr<T>& rhs ): \n")
|
||||
FILE(APPEND "${FILENAME}" " obj(rhs.get()), count(rhs.count) { if ( count!=NULL ) { ++(*count); } } \n")
|
||||
FILE(APPEND "${FILENAME}" " obj(rhs.get()), deleter(reinterpret_cast<D>(rhs.deleter)), count(rhs.count) { if ( count!=NULL ) { ++(*count); } } \n")
|
||||
FILE(APPEND "${FILENAME}" " template<class U> shared_ptr( const shared_ptr<U>& rhs ): \n")
|
||||
FILE(APPEND "${FILENAME}" " obj(rhs.get()), count(rhs.count) { if ( count!=NULL ) { ++(*count); } } \n")
|
||||
FILE(APPEND "${FILENAME}" " obj(rhs.get()), deleter(reinterpret_cast<D>(rhs.deleter)), count(rhs.count) { if ( count!=NULL ) { ++(*count); } } \n")
|
||||
FILE(APPEND "${FILENAME}" " shared_ptr& operator=( const shared_ptr<T>& rhs )\n")
|
||||
FILE(APPEND "${FILENAME}" " { obj=rhs.obj; count=rhs.count; ++(*count); return *this; } \n")
|
||||
FILE(APPEND "${FILENAME}" " { if (this==&rhs) { return *this;} reset(); obj=rhs.obj; deleter=reinterpret_cast<D>(rhs.deleter); count=rhs.count; ++(*count); return *this; } \n")
|
||||
FILE(APPEND "${FILENAME}" " ~shared_ptr( ) { reset(); }\n")
|
||||
FILE(APPEND "${FILENAME}" " void reset( T *ptr ) { reset(); obj=ptr; count=new int; (*count)=1; }\n")
|
||||
FILE(APPEND "${FILENAME}" " void reset( void ) { \n")
|
||||
FILE(APPEND "${FILENAME}" " if ( count!=NULL) { int tmp=--(*count); if ( tmp==0 ) { delete obj; delete count; } } \n")
|
||||
FILE(APPEND "${FILENAME}" " if ( count!=NULL) { int tmp=--(*count); if ( tmp==0 ) { deleter(obj); delete count; } } \n")
|
||||
FILE(APPEND "${FILENAME}" " obj=NULL; count=NULL; \n")
|
||||
FILE(APPEND "${FILENAME}" " }\n")
|
||||
FILE(APPEND "${FILENAME}" " T* get( ) const { return obj; } \n")
|
||||
|
|
|
@ -8,13 +8,13 @@
|
|||
#include "common/Utilities.h"
|
||||
|
||||
|
||||
#define GET_ARRAY_INDEX(i1,i2,i3) i1+d_N[0]*(i2+d_N[1]*i3)
|
||||
#define GET_ARRAY_INDEX(i1,i2,i3,i4) i1+d_N[0]*(i2+d_N[1]*(i3+d_N[2]*i4))
|
||||
#ifdef DEBUG
|
||||
#define CHECK_ARRAY_INDEX(i1,i2,i3) \
|
||||
if ( GET_ARRAY_INDEX(i1,i2,i3)<0 || GET_ARRAY_INDEX(i1,i2,i3)>d_length ) \
|
||||
#define CHECK_ARRAY_INDEX(i1,i2,i3,i4) \
|
||||
if ( GET_ARRAY_INDEX(i1,i2,i3,i4)<0 || GET_ARRAY_INDEX(i1,i2,i3,i4)>d_length ) \
|
||||
ERROR("Index exceeds array bounds");
|
||||
#else
|
||||
#define CHECK_ARRAY_INDEX(i1,i2,i3)
|
||||
#define CHECK_ARRAY_INDEX(i1,i2,i3,i4)
|
||||
#endif
|
||||
|
||||
|
||||
|
@ -236,42 +236,56 @@ public:
|
|||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
inline TYPE& operator()( size_t i ) { CHECK_ARRAY_INDEX(i,0,0) return d_data[i]; }
|
||||
inline TYPE& operator()( size_t i ) { CHECK_ARRAY_INDEX(i,0,0,0) return d_data[i]; }
|
||||
|
||||
/*!
|
||||
* Access the desired element
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
inline const TYPE& operator()( size_t i ) const { CHECK_ARRAY_INDEX(i,0,0) return d_data[i]; }
|
||||
inline const TYPE& operator()( size_t i ) const { CHECK_ARRAY_INDEX(i,0,0,0) return d_data[i]; }
|
||||
|
||||
/*!
|
||||
* Access the desired element
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
inline TYPE& operator()( size_t i, size_t j ) { CHECK_ARRAY_INDEX(i,j,0) return d_data[i+j*d_N[0]]; }
|
||||
inline TYPE& operator()( size_t i, size_t j ) { CHECK_ARRAY_INDEX(i,j,0,0) return d_data[i+j*d_N[0]]; }
|
||||
|
||||
/*!
|
||||
* Access the desired element
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
inline const TYPE& operator()( size_t i, size_t j ) const { CHECK_ARRAY_INDEX(i,j,0) return d_data[i+j*d_N[0]]; }
|
||||
inline const TYPE& operator()( size_t i, size_t j ) const { CHECK_ARRAY_INDEX(i,j,0,0) return d_data[i+j*d_N[0]]; }
|
||||
|
||||
/*!
|
||||
* Access the desired element
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
inline TYPE& operator()( size_t i, size_t j, size_t k ) { CHECK_ARRAY_INDEX(i,j,k) return d_data[GET_ARRAY_INDEX(i,j,k)]; }
|
||||
inline TYPE& operator()( size_t i, size_t j, size_t k ) { CHECK_ARRAY_INDEX(i,j,k,0) return d_data[GET_ARRAY_INDEX(i,j,k,0)]; }
|
||||
|
||||
/*!
|
||||
* Access the desired element
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
inline const TYPE& operator()( size_t i, size_t j, size_t k ) const { CHECK_ARRAY_INDEX(i,j,k) return d_data[GET_ARRAY_INDEX(i,j,k)]; }
|
||||
inline const TYPE& operator()( size_t i, size_t j, size_t k ) const { CHECK_ARRAY_INDEX(i,j,k,0) return d_data[GET_ARRAY_INDEX(i,j,k,0)]; }
|
||||
|
||||
/*!
|
||||
* Access the desired element
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
inline TYPE& operator()( size_t i, size_t j, size_t k, size_t m ) { CHECK_ARRAY_INDEX(i,j,k,m) return d_data[GET_ARRAY_INDEX(i,j,k,m)]; }
|
||||
|
||||
/*!
|
||||
* Access the desired element
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
inline const TYPE& operator()( size_t i, size_t j, size_t k, size_t m ) const { CHECK_ARRAY_INDEX(i,j,k,m) return d_data[GET_ARRAY_INDEX(i,j,k,m)]; }
|
||||
|
||||
|
||||
//! Check if two matricies are equal
|
||||
|
@ -312,7 +326,7 @@ public:
|
|||
|
||||
private:
|
||||
int d_ndim;
|
||||
size_t d_N[3];
|
||||
size_t d_N[4];
|
||||
size_t d_length;
|
||||
TYPE *d_data;
|
||||
std::shared_ptr<TYPE> d_ptr;
|
||||
|
|
|
@ -69,9 +69,8 @@ void Array<TYPE>::allocate( const std::vector<size_t>& N )
|
|||
d_N[0] = 0;
|
||||
d_length = 0;
|
||||
}
|
||||
if ( d_length==0 )
|
||||
d_ptr.reset();
|
||||
else
|
||||
d_ptr.reset();
|
||||
if ( d_length > 0 )
|
||||
d_ptr = std::shared_ptr<TYPE>(new TYPE[d_length],DeleteArray<TYPE>);
|
||||
d_data = d_ptr.get();
|
||||
if ( d_length>0 && d_data==NULL )
|
||||
|
@ -90,7 +89,7 @@ Array<TYPE>& Array<TYPE>::operator=( const Array& rhs )
|
|||
{
|
||||
if ( this == &rhs )
|
||||
return *this;
|
||||
this->allocate( std::vector<size_t>(rhs.d_N,rhs.d_N+rhs.d_data) );
|
||||
this->allocate( rhs.size() );
|
||||
for (size_t i=0; i<d_length; i++)
|
||||
this->d_data[i] = rhs.d_data[i];
|
||||
return *this;
|
||||
|
@ -152,15 +151,17 @@ void Array<TYPE>::resize( const std::vector<size_t>& N )
|
|||
allocate(N);
|
||||
// Copy the old values
|
||||
if ( d_length > 0 ) {
|
||||
ASSERT(sizeof(d_N)/sizeof(size_t)==3);
|
||||
ASSERT(sizeof(d_N)/sizeof(size_t)==4);
|
||||
TYPE *data1 = old_data.get();
|
||||
TYPE *data2 = d_data;
|
||||
for (size_t k=0; k<std::min(N1[2],N2[2]); k++) {
|
||||
for (size_t j=0; j<std::min(N1[1],N2[1]); j++) {
|
||||
for (size_t i=0; i<std::min(N1[0],N2[0]); i++) {
|
||||
size_t index1 = i + j*N1[0] + k*N1[0]*N1[1];
|
||||
size_t index2 = i + j*N2[0] + k*N2[0]*N2[1];
|
||||
data2[index2] = data1[index1];
|
||||
for (size_t m=0; m<std::min(N1[3],N2[3]); m++) {
|
||||
for (size_t k=0; k<std::min(N1[2],N2[2]); k++) {
|
||||
for (size_t j=0; j<std::min(N1[1],N2[1]); j++) {
|
||||
for (size_t i=0; i<std::min(N1[0],N2[0]); i++) {
|
||||
size_t index1 = i + j*N1[0] + k*N1[0]*N1[1] + m*N1[0]*N1[1]*N1[2];
|
||||
size_t index2 = i + j*N2[0] + k*N2[0]*N2[1] + m*N2[0]*N2[1]*N1[2];
|
||||
data2[index2] = data1[index1];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
@ -293,8 +294,8 @@ template<class TYPE>
|
|||
template<class TYPE2>
|
||||
void Array<TYPE>::copy( const Array<TYPE2>& array )
|
||||
{
|
||||
resize( std::vector<size_t>(array.d_N,array.d_N+array.d_ndim) );
|
||||
const TYPE2 *src = array.d_data;
|
||||
resize( array.size() );
|
||||
const TYPE2 *src = array.get();
|
||||
for (size_t i=0; i<d_length; i++)
|
||||
d_data[i] = static_cast<TYPE>(src[i]);
|
||||
}
|
||||
|
|
|
@ -0,0 +1,74 @@
|
|||
#include "common/Communication.h"
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Structure to store the rank info *
|
||||
********************************************************/
|
||||
inline int getRankForBlock( int nprocx, int nprocy, int nprocz, int i, int j, int k )
|
||||
{
|
||||
int i2 = (i+nprocx)%nprocx;
|
||||
int j2 = (j+nprocy)%nprocy;
|
||||
int k2 = (k+nprocz)%nprocz;
|
||||
return i2 + j2*nprocx + k2*nprocx*nprocy;
|
||||
}
|
||||
RankInfoStruct::RankInfoStruct()
|
||||
{
|
||||
memset(this,0,sizeof(RankInfoStruct));
|
||||
}
|
||||
RankInfoStruct::RankInfoStruct( int rank0, int nprocx, int nprocy, int nprocz )
|
||||
{
|
||||
memset(this,0,sizeof(RankInfoStruct));
|
||||
nx = nprocx;
|
||||
ny = nprocy;
|
||||
nz = nprocz;
|
||||
ix = rank0%nprocx;
|
||||
jy = (rank0/nprocx)%nprocy;
|
||||
kz = rank0/(nprocx*nprocy);
|
||||
for (int i=-1; i<=1; i++) {
|
||||
for (int j=-1; j<=1; j++) {
|
||||
for (int k=-1; k<=1; k++) {
|
||||
rank[i+1][j+1][k+1] = getRankForBlock(nprocx,nprocy,nprocz,ix+i,jy+j,kz+k);
|
||||
}
|
||||
}
|
||||
}
|
||||
ASSERT(rank[1][1][1]==rank0);
|
||||
}
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Deprecated functions *
|
||||
********************************************************/
|
||||
void InitializeRanks( const int rank, const int nprocx, const int nprocy, const int nprocz,
|
||||
int& iproc, int& jproc, int& kproc,
|
||||
int& rank_x, int& rank_y, int& rank_z,
|
||||
int& rank_X, int& rank_Y, int& rank_Z,
|
||||
int& rank_xy, int& rank_XY, int& rank_xY, int& rank_Xy,
|
||||
int& rank_xz, int& rank_XZ, int& rank_xZ, int& rank_Xz,
|
||||
int& rank_yz, int& rank_YZ, int& rank_yZ, int& rank_Yz )
|
||||
{
|
||||
const RankInfoStruct data(rank,nprocx,nprocy,nprocz);
|
||||
iproc = data.ix;
|
||||
jproc = data.jy;
|
||||
kproc = data.kz;
|
||||
rank_X = data.rank[2][1][1];
|
||||
rank_x = data.rank[0][1][1];
|
||||
rank_Y = data.rank[1][2][1];
|
||||
rank_y = data.rank[1][0][1];
|
||||
rank_Z = data.rank[1][1][2];
|
||||
rank_z = data.rank[1][1][0];
|
||||
rank_XY = data.rank[2][2][1];
|
||||
rank_xy = data.rank[0][0][1];
|
||||
rank_Xy = data.rank[2][0][1];
|
||||
rank_xY = data.rank[0][2][1];
|
||||
rank_XZ = data.rank[2][1][2];
|
||||
rank_xz = data.rank[0][1][0];
|
||||
rank_Xz = data.rank[2][1][0];
|
||||
rank_xZ = data.rank[0][1][2];
|
||||
rank_YZ = data.rank[1][2][2];
|
||||
rank_yz = data.rank[1][0][0];
|
||||
rank_Yz = data.rank[1][2][0];
|
||||
rank_yZ = data.rank[1][0][2];
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -14,6 +14,76 @@
|
|||
using namespace std;
|
||||
|
||||
|
||||
/*!
|
||||
* @brief Rank info structure
|
||||
* @details Structure used to hold the ranks for the current process and it's neighbors
|
||||
*/
|
||||
struct RankInfoStruct {
|
||||
int nx; //!< The number of processors in the x direction
|
||||
int ny; //!< The number of processors in the y direction
|
||||
int nz; //!< The number of processors in the z direction
|
||||
int ix; //!< The index of the current process in the x direction
|
||||
int jy; //!< The index of the current process in the y direction
|
||||
int kz; //!< The index of the current process in the z direction
|
||||
int rank[3][3][3]; //!< The rank for the neighbor [i][j][k]
|
||||
RankInfoStruct();
|
||||
RankInfoStruct( int rank, int nprocx, int nprocy, int nprocz );
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
* @brief Communicate halo
|
||||
* @details Fill the halo cells in an array from the neighboring processes
|
||||
*/
|
||||
template<class TYPE>
|
||||
class fillHalo
|
||||
{
|
||||
public:
|
||||
/*!
|
||||
* @brief Default constructor
|
||||
* @param[in] info Rank and neighbor rank info
|
||||
* @param[in] nx Number of local cells in the x direction
|
||||
* @param[in] ny Number of local cells in the y direction
|
||||
* @param[in] nz Number of local cells in the z direction
|
||||
* @param[in] ngx Number of ghost cells in the x direction
|
||||
* @param[in] ngy Number of ghost cells in the y direction
|
||||
* @param[in] ngz Number of ghost cells in the z direction
|
||||
* @param[in] tag Initial tag to use for the communication (we will require tag:tag+26)
|
||||
* @param[in] depth Maximum depth to support
|
||||
*/
|
||||
fillHalo( const RankInfoStruct& info, int nx, int ny, int nz,
|
||||
int ngx, int ngy, int ngz, int tag, int depth,
|
||||
bool fill_face=true, bool fill_edge=true, bool fill_corner=true );
|
||||
|
||||
//! Destructor
|
||||
~fillHalo( );
|
||||
|
||||
/*!
|
||||
* @brief Communicate the halos
|
||||
|
||||
* @param[in] array The array on which we fill the halos
|
||||
*/
|
||||
void fill( Array<TYPE>& array );
|
||||
|
||||
private:
|
||||
RankInfoStruct info;
|
||||
int nx, ny, nz, ngx, ngy, ngz, depth;
|
||||
bool fill_pattern[3][3][3];
|
||||
int tag[3][3][3];
|
||||
int N_send_recv[3][3][3];
|
||||
TYPE *mem;
|
||||
TYPE *send[3][3][3], *recv[3][3][3];
|
||||
MPI_Request send_req[3][3][3], recv_req[3][3][3];
|
||||
MPI_Comm comm;
|
||||
MPI_Datatype datatype;
|
||||
fillHalo(); // Private empty constructor
|
||||
fillHalo(const fillHalo&); // Private copy constructor
|
||||
fillHalo& operator=(const fillHalo&); // Private assignment operator
|
||||
void pack( const Array<TYPE>& array, int i, int j, int k, TYPE *buffer );
|
||||
void unpack( Array<TYPE>& array, int i, int j, int k, const TYPE *buffer );
|
||||
};
|
||||
|
||||
|
||||
//***************************************************************************************
|
||||
inline void PackMeshData(int *list, int count, double *sendbuf, double *data){
|
||||
// Fill in the phase ID values from neighboring processors
|
||||
|
@ -35,50 +105,15 @@ inline void UnpackMeshData(int *list, int count, double *recvbuf, double *data){
|
|||
}
|
||||
}
|
||||
|
||||
//***************************************************************************************
|
||||
inline int getRankForBlock( int nprocx, int nprocy, int nprocz, int i, int j, int k )
|
||||
{
|
||||
int i2 = (i+nprocx)%nprocx;
|
||||
int j2 = (j+nprocy)%nprocy;
|
||||
int k2 = (k+nprocz)%nprocz;
|
||||
return i2 + j2*nprocx + k2*nprocx*nprocy;
|
||||
}
|
||||
inline void InitializeRanks( const int rank, const int nprocx, const int nprocy, const int nprocz,
|
||||
|
||||
// Initialize the ranks (this is deprecated, see RankInfoStruct)
|
||||
void InitializeRanks( const int rank, const int nprocx, const int nprocy, const int nprocz,
|
||||
int& iproc, int& jproc, int& kproc,
|
||||
int& rank_x, int& rank_y, int& rank_z,
|
||||
int& rank_X, int& rank_Y, int& rank_Z,
|
||||
int& rank_xy, int& rank_XY, int& rank_xY, int& rank_Xy,
|
||||
int& rank_xz, int& rank_XZ, int& rank_xZ, int& rank_Xz,
|
||||
int& rank_yz, int& rank_YZ, int& rank_yZ, int& rank_Yz )
|
||||
{
|
||||
// map the rank to the block index
|
||||
iproc = rank%nprocx;
|
||||
jproc = (rank/nprocx)%nprocy;
|
||||
kproc = rank/(nprocx*nprocy);
|
||||
|
||||
// set up the neighbor ranks
|
||||
int i = iproc;
|
||||
int j = jproc;
|
||||
int k = kproc;
|
||||
rank_X = getRankForBlock(nprocx,nprocy,nprocz,i+1,j,k);
|
||||
rank_x = getRankForBlock(nprocx,nprocy,nprocz,i-1,j,k);
|
||||
rank_Y = getRankForBlock(nprocx,nprocy,nprocz,i,j+1,k);
|
||||
rank_y = getRankForBlock(nprocx,nprocy,nprocz,i,j-1,k);
|
||||
rank_Z = getRankForBlock(nprocx,nprocy,nprocz,i,j,k+1);
|
||||
rank_z = getRankForBlock(nprocx,nprocy,nprocz,i,j,k-1);
|
||||
rank_XY = getRankForBlock(nprocx,nprocy,nprocz,i+1,j+1,k);
|
||||
rank_xy = getRankForBlock(nprocx,nprocy,nprocz,i-1,j-1,k);
|
||||
rank_Xy = getRankForBlock(nprocx,nprocy,nprocz,i+1,j-1,k);
|
||||
rank_xY = getRankForBlock(nprocx,nprocy,nprocz,i-1,j+1,k);
|
||||
rank_XZ = getRankForBlock(nprocx,nprocy,nprocz,i+1,j,k+1);
|
||||
rank_xz = getRankForBlock(nprocx,nprocy,nprocz,i-1,j,k-1);
|
||||
rank_Xz = getRankForBlock(nprocx,nprocy,nprocz,i+1,j,k-1);
|
||||
rank_xZ = getRankForBlock(nprocx,nprocy,nprocz,i-1,j,k+1);
|
||||
rank_YZ = getRankForBlock(nprocx,nprocy,nprocz,i,j+1,k+1);
|
||||
rank_yz = getRankForBlock(nprocx,nprocy,nprocz,i,j-1,k-1);
|
||||
rank_Yz = getRankForBlock(nprocx,nprocy,nprocz,i,j+1,k-1);
|
||||
rank_yZ = getRankForBlock(nprocx,nprocy,nprocz,i,j-1,k+1);
|
||||
}
|
||||
int& rank_yz, int& rank_YZ, int& rank_yZ, int& rank_Yz );
|
||||
|
||||
|
||||
//***************************************************************************************
|
||||
|
@ -324,3 +359,6 @@ inline void CommunicateMeshHalo(DoubleArray &Mesh, MPI_Comm Communicator,
|
|||
|
||||
#endif
|
||||
|
||||
#include "common/Communication.hpp"
|
||||
|
||||
|
||||
|
|
202
common/Communication.hpp
Normal file
202
common/Communication.hpp
Normal file
|
@ -0,0 +1,202 @@
|
|||
#ifndef COMMUNICATION_HPP_INC
|
||||
#define COMMUNICATION_HPP_INC
|
||||
|
||||
#include "common/Communication.h"
|
||||
#include "common/MPI_Helpers.h"
|
||||
#include "ProfilerApp.h"
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Structure to store the rank info *
|
||||
********************************************************/
|
||||
template<class TYPE>
|
||||
fillHalo<TYPE>::fillHalo( const RankInfoStruct& info0, int nx0, int ny0, int nz0,
|
||||
int ngx0, int ngy0, int ngz0, int tag0, int depth0,
|
||||
bool fill_face, bool fill_edge, bool fill_corner ):
|
||||
info(info0), nx(nx0), ny(ny0), nz(nz0), ngx(ngx0), ngy(ngy0), ngz(ngz0), depth(depth0)
|
||||
{
|
||||
comm = MPI_COMM_WORLD;
|
||||
datatype = getMPItype<TYPE>();
|
||||
// Set the fill pattern
|
||||
memset(fill_pattern,0,sizeof(fill_pattern));
|
||||
if ( fill_face ) {
|
||||
fill_pattern[0][1][1] = true;
|
||||
fill_pattern[2][1][1] = true;
|
||||
fill_pattern[1][0][1] = true;
|
||||
fill_pattern[1][2][1] = true;
|
||||
fill_pattern[1][1][0] = true;
|
||||
fill_pattern[1][1][2] = true;
|
||||
}
|
||||
if ( fill_edge ) {
|
||||
fill_pattern[0][0][1] = true;
|
||||
fill_pattern[0][2][1] = true;
|
||||
fill_pattern[2][0][1] = true;
|
||||
fill_pattern[2][2][1] = true;
|
||||
fill_pattern[0][1][0] = true;
|
||||
fill_pattern[0][1][2] = true;
|
||||
fill_pattern[2][1][0] = true;
|
||||
fill_pattern[2][1][2] = true;
|
||||
fill_pattern[1][0][0] = true;
|
||||
fill_pattern[1][0][2] = true;
|
||||
fill_pattern[1][2][0] = true;
|
||||
fill_pattern[1][2][2] = true;
|
||||
}
|
||||
if ( fill_corner ) {
|
||||
fill_pattern[0][0][0] = true;
|
||||
fill_pattern[0][0][2] = true;
|
||||
fill_pattern[0][2][0] = true;
|
||||
fill_pattern[0][2][2] = true;
|
||||
fill_pattern[2][0][0] = true;
|
||||
fill_pattern[2][0][2] = true;
|
||||
fill_pattern[2][2][0] = true;
|
||||
fill_pattern[2][2][2] = true;
|
||||
}
|
||||
// Determine the number of elements for each send/recv
|
||||
for (int i=0; i<3; i++) {
|
||||
int ni = (i-1)==0 ? nx:ngx;
|
||||
for (int j=0; j<3; j++) {
|
||||
int nj = (j-1)==0 ? ny:ngy;
|
||||
for (int k=0; k<3; k++) {
|
||||
int nk = (k-1)==0 ? nz:ngz;
|
||||
if ( fill_pattern[i][j][k] )
|
||||
N_send_recv[i][j][k] = ni*nj*nk;
|
||||
else
|
||||
N_send_recv[i][j][k] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Create send/recv buffers
|
||||
size_t N_mem=0;
|
||||
for (int i=0; i<3; i++) {
|
||||
for (int j=0; j<3; j++) {
|
||||
for (int k=0; k<3; k++)
|
||||
N_mem += N_send_recv[i][j][k];
|
||||
}
|
||||
}
|
||||
mem = new TYPE[2*depth*N_mem];
|
||||
size_t index = 0;
|
||||
for (int i=0; i<3; i++) {
|
||||
for (int j=0; j<3; j++) {
|
||||
for (int k=0; k<3; k++) {
|
||||
send[i][j][k] = &mem[index];
|
||||
index += depth*N_send_recv[i][j][k];
|
||||
recv[i][j][k] = &mem[index];
|
||||
index += depth*N_send_recv[i][j][k];
|
||||
}
|
||||
}
|
||||
}
|
||||
// Create the tags
|
||||
for (int i=0; i<3; i++) {
|
||||
for (int j=0; j<3; j++) {
|
||||
for (int k=0; k<3; k++) {
|
||||
tag[i][j][k] = tag0 + i + j*3 + k*9;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
template<class TYPE>
|
||||
fillHalo<TYPE>::~fillHalo( )
|
||||
{
|
||||
delete [] mem;
|
||||
}
|
||||
template<class TYPE>
|
||||
void fillHalo<TYPE>::fill( Array<TYPE>& data )
|
||||
{
|
||||
PROFILE_START("fillHalo::fill",1);
|
||||
int depth2 = data.size(3);
|
||||
ASSERT((int)data.size(0)==nx+2*ngx);
|
||||
ASSERT((int)data.size(1)==ny+2*ngy);
|
||||
ASSERT((int)data.size(2)==nz+2*ngz);
|
||||
ASSERT(depth2<=depth);
|
||||
ASSERT(data.ndim()==3||data.ndim()==4);
|
||||
// Start the recieves
|
||||
for (int i=0; i<3; i++) {
|
||||
for (int j=0; j<3; j++) {
|
||||
for (int k=0; k<3; k++) {
|
||||
if ( !fill_pattern[i][j][k] )
|
||||
continue;
|
||||
MPI_Irecv( recv[i][j][k], depth2*N_send_recv[i][j][k], datatype,
|
||||
info.rank[i][j][k], tag[2-i][2-j][2-k], comm, &recv_req[i][j][k] );
|
||||
}
|
||||
}
|
||||
}
|
||||
// Pack the src data and start the sends
|
||||
for (int i=0; i<3; i++) {
|
||||
for (int j=0; j<3; j++) {
|
||||
for (int k=0; k<3; k++) {
|
||||
if ( !fill_pattern[i][j][k] )
|
||||
continue;
|
||||
pack( data, i-1, j-1, k-1, send[i][j][k] );
|
||||
MPI_Isend( send[i][j][k], depth2*N_send_recv[i][j][k], datatype,
|
||||
info.rank[i][j][k], tag[i][j][k], comm, &send_req[i][j][k] );
|
||||
}
|
||||
}
|
||||
}
|
||||
// Recv the dst data and unpack
|
||||
MPI_Status status;
|
||||
for (int i=0; i<3; i++) {
|
||||
for (int j=0; j<3; j++) {
|
||||
for (int k=0; k<3; k++) {
|
||||
if ( !fill_pattern[i][j][k] )
|
||||
continue;
|
||||
MPI_Wait(&recv_req[i][j][k],&status);
|
||||
unpack( data, i-1, j-1, k-1, recv[i][j][k] );
|
||||
}
|
||||
}
|
||||
}
|
||||
// Wait until all sends have completed
|
||||
for (int i=0; i<3; i++) {
|
||||
for (int j=0; j<3; j++) {
|
||||
for (int k=0; k<3; k++) {
|
||||
if ( !fill_pattern[i][j][k] )
|
||||
continue;
|
||||
MPI_Wait(&send_req[i][j][k],&status);
|
||||
}
|
||||
}
|
||||
}
|
||||
PROFILE_STOP("fillHalo::fill",1);
|
||||
}
|
||||
template<class TYPE>
|
||||
void fillHalo<TYPE>::pack( const Array<TYPE>& data, int i0, int j0, int k0, TYPE *buffer )
|
||||
{
|
||||
int depth2 = data.size(3);
|
||||
int ni = i0==0 ? nx:ngx;
|
||||
int nj = j0==0 ? ny:ngy;
|
||||
int nk = k0==0 ? nz:ngz;
|
||||
int is = i0==0 ? ngx:((i0==-1)?ngx:nx);
|
||||
int js = j0==0 ? ngy:((j0==-1)?ngy:ny);
|
||||
int ks = k0==0 ? ngz:((k0==-1)?ngz:nz);
|
||||
for (int d=0; d<depth2; d++) {
|
||||
for (int k=0; k<nk; k++) {
|
||||
for (int j=0; j<nj; j++) {
|
||||
for (int i=0; i<ni; i++) {
|
||||
buffer[i+j*ni+k*ni*nj+d*ni*nj*nk] = data(i+is,j+js,k+ks,d);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
template<class TYPE>
|
||||
void fillHalo<TYPE>::unpack( Array<TYPE>& data, int i0, int j0, int k0, const TYPE *buffer )
|
||||
{
|
||||
int depth2 = data.size(3);
|
||||
int ni = i0==0 ? nx:ngx;
|
||||
int nj = j0==0 ? ny:ngy;
|
||||
int nk = k0==0 ? nz:ngz;
|
||||
int is = i0==0 ? ngx:((i0==-1)?0:nx+ngx);
|
||||
int js = j0==0 ? ngy:((j0==-1)?0:ny+ngy);
|
||||
int ks = k0==0 ? ngz:((k0==-1)?0:nz+ngz);
|
||||
for (int d=0; d<depth2; d++) {
|
||||
for (int k=0; k<nk; k++) {
|
||||
for (int j=0; j<nj; j++) {
|
||||
for (int i=0; i<ni; i++) {
|
||||
data(i+is,j+js,k+ks,d) = buffer[i+j*ni+k*ni*nj+d*ni*nj*nk];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
|
@ -2,6 +2,34 @@
|
|||
#include "common/Utilities.h"
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Return the MPI data type *
|
||||
********************************************************/
|
||||
template<> MPI_Datatype getMPItype<char>() {
|
||||
return MPI_CHAR;
|
||||
}
|
||||
template<> MPI_Datatype getMPItype<unsigned char>() {
|
||||
return MPI_UNSIGNED_CHAR;
|
||||
}
|
||||
template<> MPI_Datatype getMPItype<int>() {
|
||||
return MPI_INT;
|
||||
}
|
||||
template<> MPI_Datatype getMPItype<long>() {
|
||||
return MPI_LONG;
|
||||
}
|
||||
template<> MPI_Datatype getMPItype<unsigned long>() {
|
||||
return MPI_UNSIGNED_LONG;
|
||||
}
|
||||
template<> MPI_Datatype getMPItype<long long>() {
|
||||
return MPI_LONG_LONG;
|
||||
}
|
||||
template<> MPI_Datatype getMPItype<float>() {
|
||||
return MPI_FLOAT;
|
||||
}
|
||||
template<> MPI_Datatype getMPItype<double>() {
|
||||
return MPI_DOUBLE;
|
||||
}
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Concrete implimentations for packing/unpacking *
|
||||
|
|
|
@ -57,6 +57,10 @@ inline int MPI_WORLD_RANK( ) {
|
|||
return rank;
|
||||
}
|
||||
|
||||
//! Return the appropriate MPI datatype for a class
|
||||
template<class TYPE>
|
||||
MPI_Datatype getMPItype();
|
||||
|
||||
|
||||
//! Template function to return the buffer size required to pack a class
|
||||
template<class TYPE>
|
||||
|
|
|
@ -10,7 +10,6 @@
|
|||
|
||||
|
||||
|
||||
|
||||
/********************************************************
|
||||
* Default instantiations for std::vector *
|
||||
********************************************************/
|
||||
|
|
132
common/pmmc.h
132
common/pmmc.h
|
@ -370,138 +370,6 @@ public:
|
|||
}
|
||||
|
||||
};
|
||||
//--------------------------------------------------------------------------------------------------------
|
||||
inline int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator,
|
||||
DoubleArray &F, DoubleArray &S, double vf, double vs, int startx, int starty,
|
||||
int startz, IntArray &temp)
|
||||
{
|
||||
// Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob
|
||||
// F>vf => oil phase S>vs => in porespace
|
||||
// update the list of blobs, indicator mesh
|
||||
int m = F.size(0); // maxima for the meshes
|
||||
int n = F.size(1);
|
||||
int o = F.size(2);
|
||||
|
||||
int cubes_in_blob=0;
|
||||
int nrecent = 1; // number of nodes added at most recent sweep
|
||||
temp(0,0) = startx; // Set the initial point as a "seed" for the sweeps
|
||||
temp(1,0) = starty;
|
||||
temp(2,0) = startz;
|
||||
int ntotal = 1; // total number of nodes in blob
|
||||
indicator(startx,starty,startz) = nblobs;
|
||||
|
||||
int p,s,x,y,z,start,finish,nodx,nody,nodz;
|
||||
int imin=startx,imax=startx,jmin=starty,jmax=starty; // initialize maxima / minima
|
||||
int kmin=startz,kmax=startz;
|
||||
int d[26][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},
|
||||
{1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0},{1,0,1},{-1,0,1},
|
||||
{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-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}}; // directions to neighbors
|
||||
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
|
||||
bool status = 1; // status == true => continue to look for points
|
||||
while (status == 1){
|
||||
start = ntotal - nrecent;
|
||||
finish = ntotal;
|
||||
nrecent = 0; // set recent points back to zero for next sweep through
|
||||
for (s=start;s<finish;s++){
|
||||
// Loop over recent points; look for new points
|
||||
x = temp(0,s);
|
||||
y = temp(1,s);
|
||||
z = temp(2,s);
|
||||
// Looop over the directions
|
||||
for (p=0;p<26;p++){
|
||||
nodx=x+d[p][0];
|
||||
|
||||
if (nodx < 0 ){ nodx = m-1; } // Periodic BC for x
|
||||
if (nodx > m-1 ){ nodx = 0; }
|
||||
nody=y+d[p][1];
|
||||
|
||||
if (nody < 0 ){ nody = n-1; } // Periodic BC for y
|
||||
if (nody > n-1 ){ nody = 0; }
|
||||
|
||||
nodz=z+d[p][2];
|
||||
if (nodz < 0 ){ nodz = o-1; } // Periodic BC for z
|
||||
if (nodz > o-1 ){ nodz = 0; }
|
||||
|
||||
if ( F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs
|
||||
&& indicator(nodx,nody,nodz) == -1 ){
|
||||
// Node is a part of the blob - add it to the list
|
||||
temp(0,ntotal) = nodx;
|
||||
temp(1,ntotal) = nody;
|
||||
temp(2,ntotal) = nodz;
|
||||
ntotal++;
|
||||
nrecent++;
|
||||
// Update the indicator map
|
||||
indicator(nodx,nody,nodz) = nblobs;
|
||||
// Update the min / max for the cube loop
|
||||
if ( nodx < imin ){ imin = nodx; }
|
||||
if ( nodx > imax ){ imax = nodx; }
|
||||
if ( nody < jmin ){ jmin = nody; }
|
||||
if ( nody > jmax ){ jmax = nody; }
|
||||
if ( nodz < kmin ){ kmin = nodz; }
|
||||
if ( nodz > kmax ){ kmax = nodz; }
|
||||
}
|
||||
else if (F(nodx,nody,nodz) > vf && S(nodx,nody,nodz) > vs
|
||||
&& indicator(nodx,nody,nodz) > -1 && indicator(nodx,nody,nodz) != nblobs){
|
||||
// Some kind of error in algorithm
|
||||
printf("Error in blob search algorithm!");
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
if ( nrecent == 0){
|
||||
status = 0;
|
||||
}
|
||||
}
|
||||
// Use points in temporary storage array to add cubes to the list of blobs
|
||||
if ( imin > 0) { imin = imin-1; }
|
||||
// if ( imax < m-1) { imax = imax+1; }
|
||||
if ( jmin > 0) { jmin = jmin-1; }
|
||||
// if ( jmax < n-1) { jmax = jmax+1; }
|
||||
if ( kmin > 0) { kmin = kmin-1; }
|
||||
// if ( kmax < o-1) { kmax = kmax+1; }
|
||||
int i,j,k;
|
||||
bool add;
|
||||
for (k=kmin;k<kmax;k++){
|
||||
for (j=jmin;j<jmax;j++){
|
||||
for (i=imin;i<imax;i++){
|
||||
// If cube(i,j,k) has any nodes in blob, add it to the list
|
||||
// Loop over cube edges
|
||||
add = 0;
|
||||
for (p=0;p<8;p++){
|
||||
nodx = i+cube[p][0];
|
||||
nody = j+cube[p][1];
|
||||
nodz = k+cube[p][2];
|
||||
if ( indicator(nodx,nody,nodz) == nblobs ){
|
||||
// Cube corner is in this blob
|
||||
add = 1;
|
||||
}
|
||||
}
|
||||
if (add == 1){
|
||||
// Add cube to the list
|
||||
blobs(0,ncubes) = i;
|
||||
blobs(1,ncubes) = j;
|
||||
blobs(2,ncubes) = k;
|
||||
ncubes++;
|
||||
cubes_in_blob++;
|
||||
// Loop again to check for overlap
|
||||
for (p=0;p<8;p++){
|
||||
nodx = i+cube[p][0];
|
||||
nody = j+cube[p][1];
|
||||
nodz = k+cube[p][2];
|
||||
if (indicator(nodx,nody,nodz) > -1 && indicator(nodx,nody,nodz) != nblobs){
|
||||
printf("Overlapping cube!");
|
||||
cout << i << ", " << j << ", " << k << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return cubes_in_blob;
|
||||
}
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/*
|
||||
inline DoubleArray SOLVE( DoubleArray &A, DoubleArray &b)
|
||||
|
|
|
@ -5,7 +5,8 @@
|
|||
|
||||
#include <iostream>
|
||||
#include <math.h>
|
||||
#include "pmmc.h"
|
||||
#include "common/pmmc.h"
|
||||
#include "analysis/analysis.h"
|
||||
//#include "Domain.h"
|
||||
|
||||
#define NUM_AVERAGES 30
|
||||
|
|
|
@ -5,436 +5,336 @@
|
|||
|
||||
#include <iostream>
|
||||
#include <math.h>
|
||||
#include "pmmc.h"
|
||||
#include "common/pmmc.h"
|
||||
#include "analysis/analysis.h"
|
||||
|
||||
//#include "Domain.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
inline void ReadCheckpoint(char *FILENAME, double *cDen, double *cDistEven, double *cDistOdd, int N)
|
||||
{
|
||||
int q,n;
|
||||
double value;
|
||||
ifstream File(FILENAME,ios::binary);
|
||||
for (n=0; n<N; n++){
|
||||
// Write the two density values
|
||||
File.read((char*) &value, sizeof(value));
|
||||
cDen[n] = value;
|
||||
// if (n== 66276) printf("Density a = %f \n",value);
|
||||
File.read((char*) &value, sizeof(value));
|
||||
cDen[N+n] = value;
|
||||
// if (n== 66276) printf("Density b = %f \n",value);
|
||||
// Read the even distributions
|
||||
for (q=0; q<10; q++){
|
||||
File.read((char*) &value, sizeof(value));
|
||||
cDistEven[q*N+n] = value;
|
||||
// if (n== 66276) printf("dist even %i = %f \n",q,value);
|
||||
}
|
||||
// Read the odd distributions
|
||||
for (q=0; q<9; q++){
|
||||
File.read((char*) &value, sizeof(value));
|
||||
cDistOdd[q*N+n] = value;
|
||||
// if (n== 66276) printf("dist even %i = %f \n",q,value);
|
||||
}
|
||||
}
|
||||
File.close();
|
||||
int q,n;
|
||||
double value;
|
||||
ifstream File(FILENAME,ios::binary);
|
||||
for (n=0; n<N; n++){
|
||||
// Write the two density values
|
||||
File.read((char*) &value, sizeof(value));
|
||||
cDen[n] = value;
|
||||
// if (n== 66276) printf("Density a = %f \n",value);
|
||||
File.read((char*) &value, sizeof(value));
|
||||
cDen[N+n] = value;
|
||||
// if (n== 66276) printf("Density b = %f \n",value);
|
||||
// Read the even distributions
|
||||
for (q=0; q<10; q++){
|
||||
File.read((char*) &value, sizeof(value));
|
||||
cDistEven[q*N+n] = value;
|
||||
// if (n== 66276) printf("dist even %i = %f \n",q,value);
|
||||
}
|
||||
// Read the odd distributions
|
||||
for (q=0; q<9; q++){
|
||||
File.read((char*) &value, sizeof(value));
|
||||
cDistOdd[q*N+n] = value;
|
||||
// if (n== 66276) printf("dist even %i = %f \n",q,value);
|
||||
}
|
||||
}
|
||||
File.close();
|
||||
}
|
||||
|
||||
inline void ReadBinaryFile(char *FILENAME, double *Data, int N)
|
||||
{
|
||||
int n;
|
||||
double value;
|
||||
ifstream File(FILENAME,ios::binary);
|
||||
for (n=0; n<N; n++){
|
||||
// Write the two density values
|
||||
File.read((char*) &value, sizeof(value));
|
||||
Data[n] = value;
|
||||
int n;
|
||||
double value;
|
||||
ifstream File(FILENAME,ios::binary);
|
||||
for (n=0; n<N; n++){
|
||||
// Write the two density values
|
||||
File.read((char*) &value, sizeof(value));
|
||||
Data[n] = value;
|
||||
|
||||
}
|
||||
File.close();
|
||||
}
|
||||
File.close();
|
||||
}
|
||||
|
||||
inline void SetPeriodicBC(DoubleArray &Scalar, int nx, int ny, int nz){
|
||||
|
||||
int i,j,k,in,jn,kn;
|
||||
for (k=0; k<nz; k++){
|
||||
for (j=0; j<ny; j++){
|
||||
for (i=0; i<nx; i++){
|
||||
in = i; jn=j; kn=k;
|
||||
if (i==0) in = nx-2 ;
|
||||
else if (i==nx-1) in = 0;
|
||||
if (j==0) jn = ny-2;
|
||||
else if (j==ny-1) jn = 0;
|
||||
if (k==0) kn = nz-2;
|
||||
else if (k==nz-1) kn = 0;
|
||||
Scalar(i,j,k) = Scalar(in,jn,kn);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int i,j,k,in,jn,kn;
|
||||
for (k=0; k<nz; k++){
|
||||
for (j=0; j<ny; j++){
|
||||
for (i=0; i<nx; i++){
|
||||
in = i; jn=j; kn=k;
|
||||
if (i==0) in = nx-2 ;
|
||||
else if (i==nx-1) in = 0;
|
||||
if (j==0) jn = ny-2;
|
||||
else if (j==ny-1) jn = 0;
|
||||
if (k==0) kn = nz-2;
|
||||
else if (k==nz-1) kn = 0;
|
||||
Scalar(i,j,k) = Scalar(in,jn,kn);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
inline void ReadFromRank(char *FILENAME, DoubleArray &Phase, int nx, int ny, int nz, int iproc, int
|
||||
jproc, int kproc)
|
||||
jproc, int kproc)
|
||||
{
|
||||
int i,j,k,q,n,N;
|
||||
int iglobal,jglobal,kglobal;
|
||||
double value;
|
||||
double denA,denB;
|
||||
int i,j,k,q,n,N;
|
||||
int iglobal,jglobal,kglobal;
|
||||
double value;
|
||||
double denA,denB;
|
||||
|
||||
N = nx*ny*nz;
|
||||
|
||||
double *Den;
|
||||
|
||||
Den = new double[2*N];
|
||||
N = nx*ny*nz;
|
||||
|
||||
double *Den;
|
||||
|
||||
Den = new double[2*N];
|
||||
|
||||
ifstream File(FILENAME,ios::binary);
|
||||
for (n=0; n<N; n++){
|
||||
// Write the two density values
|
||||
File.read((char*) &value, sizeof(value));
|
||||
Den[2*n] = value;
|
||||
// if (n== 66276) printf("Density a = %f \n",value);
|
||||
File.read((char*) &value, sizeof(value));
|
||||
Den[2*n+1] = value;
|
||||
ifstream File(FILENAME,ios::binary);
|
||||
for (n=0; n<N; n++){
|
||||
// Write the two density values
|
||||
File.read((char*) &value, sizeof(value));
|
||||
Den[2*n] = value;
|
||||
// if (n== 66276) printf("Density a = %f \n",value);
|
||||
File.read((char*) &value, sizeof(value));
|
||||
Den[2*n+1] = value;
|
||||
|
||||
// if (n== 66276) printf("Density b = %f \n",value);
|
||||
// Read the even distributions
|
||||
for (q=0; q<10; q++){
|
||||
File.read((char*) &value, sizeof(value));
|
||||
}
|
||||
// Read the odd distributions
|
||||
for (q=0; q<9; q++){
|
||||
File.read((char*) &value, sizeof(value));
|
||||
}
|
||||
}
|
||||
File.close();
|
||||
|
||||
// Compute the phase field
|
||||
for (k=1; k<nz-1; k++){
|
||||
for (j=1; j<ny-1; j++){
|
||||
for (i=1; i<nz-1; i++){
|
||||
//........................................................................
|
||||
n = k*nx*ny+j*nx+i;
|
||||
//........................................................................
|
||||
denA = Den[n];
|
||||
denB = Den[N+n];
|
||||
//........................................................................
|
||||
// save values in global arrays
|
||||
//........................................................................
|
||||
iglobal = iproc*(nx-2)+i;
|
||||
jglobal = jproc*(ny-2)+j;
|
||||
kglobal = kproc*(nz-2)+k;
|
||||
//........................................................................
|
||||
Phase(iglobal,jglobal,kglobal) = (denA-denB)/(denA+denB);
|
||||
//........................................................................
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
delete Den;
|
||||
// if (n== 66276) printf("Density b = %f \n",value);
|
||||
// Read the even distributions
|
||||
for (q=0; q<10; q++){
|
||||
File.read((char*) &value, sizeof(value));
|
||||
}
|
||||
// Read the odd distributions
|
||||
for (q=0; q<9; q++){
|
||||
File.read((char*) &value, sizeof(value));
|
||||
}
|
||||
}
|
||||
File.close();
|
||||
|
||||
// Compute the phase field
|
||||
for (k=1; k<nz-1; k++){
|
||||
for (j=1; j<ny-1; j++){
|
||||
for (i=1; i<nz-1; i++){
|
||||
//........................................................................
|
||||
n = k*nx*ny+j*nx+i;
|
||||
//........................................................................
|
||||
denA = Den[n];
|
||||
denB = Den[N+n];
|
||||
//........................................................................
|
||||
// save values in global arrays
|
||||
//........................................................................
|
||||
iglobal = iproc*(nx-2)+i;
|
||||
jglobal = jproc*(ny-2)+j;
|
||||
kglobal = kproc*(nz-2)+k;
|
||||
//........................................................................
|
||||
Phase(iglobal,jglobal,kglobal) = (denA-denB)/(denA+denB);
|
||||
//........................................................................
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
delete Den;
|
||||
}
|
||||
|
||||
|
||||
void readRankData( int proc, int nx, int ny, int nz, DoubleArray& Phase, DoubleArray& SignDist )
|
||||
{
|
||||
Phase.resize(nx,ny,nz);
|
||||
SignDist.resize(nx,ny,nz);
|
||||
char file1[40], file2[40];
|
||||
sprintf(file1,"SignDist.%05d",proc);
|
||||
sprintf(file2,"Phase.%05d",proc);
|
||||
ReadBinaryFile(file1, Phase.get(), nx*ny*nz);
|
||||
ReadBinaryFile(file2, SignDist.get(), nx*ny*nz);
|
||||
}
|
||||
|
||||
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
printf("-----------------------------------------------------------\n");
|
||||
printf("Labeling Blobs from Two-Phase Lattice Boltzmann Simulation \n");
|
||||
printf("-----------------------------------------------------------\n");
|
||||
// Initialize MPI
|
||||
MPI_Init(&argc,&argv);
|
||||
|
||||
//.......................................................................
|
||||
int nprocx,nprocy,nprocz,nprocs;
|
||||
int Nx, Ny, Nz;
|
||||
int nx,ny,nz;
|
||||
int nspheres;
|
||||
double Lx,Ly,Lz;
|
||||
//.......................................................................
|
||||
int i,j,k,n,p,idx;
|
||||
int iproc,jproc,kproc;
|
||||
//.......................................................................
|
||||
// Reading the domain information file
|
||||
//.......................................................................
|
||||
ifstream domain("Domain.in");
|
||||
domain >> nprocx;
|
||||
domain >> nprocy;
|
||||
domain >> nprocz;
|
||||
domain >> nx;
|
||||
domain >> ny;
|
||||
domain >> nz;
|
||||
domain >> nspheres;
|
||||
domain >> Lx;
|
||||
domain >> Ly;
|
||||
domain >> Lz;
|
||||
//.......................................................................
|
||||
printf("-----------------------------------------------------------\n");
|
||||
printf("Labeling Blobs from Two-Phase Lattice Boltzmann Simulation \n");
|
||||
printf("-----------------------------------------------------------\n");
|
||||
|
||||
nx+=2;
|
||||
ny+=2;
|
||||
nz+=2;
|
||||
|
||||
nprocs = nprocx*nprocy*nprocz;
|
||||
printf("Number of MPI ranks: %i \n", nprocs);
|
||||
Nx = (nx-2)*nprocx+2;
|
||||
Ny = (ny-2)*nprocy+2;
|
||||
Nz = (nz-2)*nprocz+2;
|
||||
printf("Full domain size: %i x %i x %i \n", Nx,Ny,Nz);
|
||||
|
||||
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
|
||||
//.......................................................................
|
||||
int nprocx,nprocy,nprocz,nprocs;
|
||||
int Nx, Ny, Nz;
|
||||
int nx,ny,nz;
|
||||
int nspheres;
|
||||
double Lx,Ly,Lz;
|
||||
|
||||
DoubleArray Phase(Nx,Ny,Nz);
|
||||
DoubleArray SignDist(Nx,Ny,Nz);
|
||||
|
||||
// Filenames used
|
||||
char LocalRankString[8];
|
||||
char LocalRankFilename[40];
|
||||
char LocalRestartFile[40];
|
||||
char BaseFilename[20];
|
||||
char tmpstr[10];
|
||||
//.......................................................................
|
||||
// Reading the domain information file
|
||||
//.......................................................................
|
||||
ifstream domain("Domain.in");
|
||||
domain >> nprocx;
|
||||
domain >> nprocy;
|
||||
domain >> nprocz;
|
||||
domain >> nx;
|
||||
domain >> ny;
|
||||
domain >> nz;
|
||||
domain >> nspheres;
|
||||
domain >> Lx;
|
||||
domain >> Ly;
|
||||
domain >> Lz;
|
||||
//.......................................................................
|
||||
|
||||
int proc,iglobal,kglobal,jglobal;
|
||||
nx+=2;
|
||||
ny+=2;
|
||||
nz+=2;
|
||||
|
||||
nprocs = nprocx*nprocy*nprocz;
|
||||
printf("Number of MPI ranks: %i \n", nprocs);
|
||||
Nx = (nx-2)*nprocx+2;
|
||||
Ny = (ny-2)*nprocy+2;
|
||||
Nz = (nz-2)*nprocz+2;
|
||||
printf("Full domain size: %i x %i x %i \n", Nx,Ny,Nz);
|
||||
|
||||
double * Temp;
|
||||
Temp = new double[nx*ny*nz];
|
||||
|
||||
for (k=0; k<Nz; k++){
|
||||
for (j=0; j<Ny; j++){
|
||||
for (i=0; i<Nx; i++){
|
||||
SignDist(i,j,k) = -100.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// read the files and populate main arrays
|
||||
for ( kproc=0; kproc<nprocz; kproc++){
|
||||
for ( jproc=0; jproc<nprocy; jproc++){
|
||||
for ( iproc=0; iproc<nprocx; iproc++){
|
||||
|
||||
proc = kproc*nprocx*nprocy + jproc*nprocx + iproc;
|
||||
DoubleArray Phase(Nx,Ny,Nz);
|
||||
DoubleArray SignDist(Nx,Ny,Nz);
|
||||
Phase.fill(0);
|
||||
SignDist.fill(0);
|
||||
|
||||
// Filenames used
|
||||
char LocalRankString[8];
|
||||
char LocalRankFilename[40];
|
||||
char LocalRestartFile[40];
|
||||
char BaseFilename[20];
|
||||
char tmpstr[10];
|
||||
|
||||
for (int k=0; k<Nz; k++){
|
||||
for (int j=0; j<Ny; j++){
|
||||
for (int i=0; i<Nx; i++){
|
||||
SignDist(i,j,k) = -100.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// read the files and populate main arrays
|
||||
for (int kproc=0; kproc<nprocz; kproc++){
|
||||
for (int jproc=0; jproc<nprocy; jproc++){
|
||||
for (int iproc=0; iproc<nprocx; iproc++){
|
||||
|
||||
sprintf(LocalRankString,"%05d",proc);
|
||||
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
|
||||
ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);
|
||||
for (k=1; k<nz-1; k++){
|
||||
for (j=1; j<ny-1; j++){
|
||||
for (i=1; i<nz-1; i++){
|
||||
int proc = kproc*nprocx*nprocy + jproc*nprocx + iproc;
|
||||
DoubleArray PhaseTmp;
|
||||
DoubleArray SignDistTmp;
|
||||
readRankData( proc, nx, ny, nz, PhaseTmp, SignDistTmp );
|
||||
|
||||
//........................................................................
|
||||
n = k*nx*ny+j*nx+i;
|
||||
//........................................................................
|
||||
iglobal = iproc*(nx-2)+i;
|
||||
jglobal = jproc*(ny-2)+j;
|
||||
kglobal = kproc*(nz-2)+k;
|
||||
//........................................................................
|
||||
SignDist(iglobal,jglobal,kglobal) = Temp[n];
|
||||
//........................................................................
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
|
||||
ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);
|
||||
for (k=1; k<nz-1; k++){
|
||||
for (j=1; j<ny-1; j++){
|
||||
for (i=1; i<nx-1; i++){
|
||||
for (int k=1; k<nz-1; k++){
|
||||
for (int j=1; j<ny-1; j++){
|
||||
for (int i=1; i<nz-1; i++){
|
||||
int iglobal = iproc*(nx-2)+i;
|
||||
int jglobal = jproc*(ny-2)+j;
|
||||
int kglobal = kproc*(nz-2)+k;
|
||||
SignDist(iglobal,jglobal,kglobal) = SignDistTmp(i,j,k);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (int k=1; k<nz-1; k++){
|
||||
for (int j=1; j<ny-1; j++){
|
||||
for (int i=1; i<nx-1; i++){
|
||||
int iglobal = iproc*(nx-2)+i;
|
||||
int jglobal = jproc*(ny-2)+j;
|
||||
int kglobal = kproc*(nz-2)+k;
|
||||
Phase(iglobal,jglobal,kglobal) = PhaseTmp(i,j,k);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
printf("Read %i ranks of Phase, SignDist \n",nprocs);
|
||||
|
||||
//........................................................................
|
||||
n = k*nx*ny+j*nx+i;
|
||||
//........................................................................
|
||||
iglobal = iproc*(nx-2)+i;
|
||||
jglobal = jproc*(ny-2)+j;
|
||||
kglobal = kproc*(nz-2)+k;
|
||||
//........................................................................
|
||||
Phase(iglobal,jglobal,kglobal) = Temp[n];
|
||||
//........................................................................
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
printf("Read %i ranks of Phase, SignDist \n",nprocs);
|
||||
|
||||
delete [] Temp;
|
||||
|
||||
IntArray GlobalBlobID(Nx,Ny,Nz);
|
||||
SetPeriodicBC(SignDist, Nx, Ny, Nz);
|
||||
SetPeriodicBC(Phase, Nx, Ny, Nz);
|
||||
|
||||
//FILE *PHASE;
|
||||
//PHASE = fopen("Phase.dat","wb");
|
||||
//fwrite(Phase.data,8,Nx*Ny*Nz,PHASE);
|
||||
//fclose(PHASE);
|
||||
|
||||
|
||||
// Compute the porosity
|
||||
double porosity=0.0;
|
||||
for (int k=0; k<Nz; k++){
|
||||
for (int j=0; j<Ny; j++){
|
||||
for (int i=0; i<Nx; i++){
|
||||
if (SignDist(i,j,k) > 0.0){
|
||||
porosity += 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
//int N=int(porosity*1.25);
|
||||
//porosity /= (Nx*Ny*Nz*1.0);
|
||||
//printf("Media porosity is %f \n",porosity);
|
||||
|
||||
SetPeriodicBC(SignDist, Nx, Ny, Nz);
|
||||
SetPeriodicBC(Phase, Nx, Ny, Nz);
|
||||
|
||||
// FILE *PHASE;
|
||||
//PHASE = fopen("Phase.dat","wb");
|
||||
//fwrite(Phase.data,8,Nx*Ny*Nz,PHASE);
|
||||
//fclose(PHASE);
|
||||
|
||||
// Initialize the local blob ID
|
||||
// Initializing the blob ID
|
||||
for (k=0; k<Nz; k++){
|
||||
for (j=0; j<Ny; j++){
|
||||
for (i=0; i<Nx; i++){
|
||||
if (SignDist(i,j,k) < 0.0){
|
||||
// Solid phase
|
||||
GlobalBlobID(i,j,k) = -2;
|
||||
}
|
||||
else{
|
||||
GlobalBlobID(i,j,k) = -1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Compute the porosity
|
||||
double porosity=0.0;
|
||||
for (k=0; k<Nz; k++){
|
||||
for (j=0; j<Ny; j++){
|
||||
for (i=0; i<Nx; i++){
|
||||
if (SignDist(i,j,k) > 0.0){
|
||||
porosity += 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
int N=int(porosity*1.25);
|
||||
porosity /= (Nx*Ny*Nz*1.0);
|
||||
printf("Media porosity is %f \n",porosity);
|
||||
|
||||
/* ****************************************************************
|
||||
IDENTIFY ALL BLOBS: F > vF, S > vS
|
||||
****************************************************************** */
|
||||
// Find blob domains, number of blobs
|
||||
int nblobs = 0; // number of blobs
|
||||
int ncubes = 0; // total number of nodes in any blob
|
||||
IntArray blobs(3,N); // store indices for blobs (cubes)
|
||||
IntArray temp(3,N); // temporary storage array
|
||||
IntArray b(N); // number of nodes in each blob
|
||||
/* ****************************************************************
|
||||
IDENTIFY ALL BLOBS: F > vF, S > vS
|
||||
****************************************************************** */
|
||||
// Find blob domains, number of blobs
|
||||
double vF=0.0;
|
||||
double vS=0.0;
|
||||
printf("Execute blob identification algorithm... \n");
|
||||
IntArray GlobalBlobID;
|
||||
int nblobs = ComputeLocalBlobIDs( Phase, SignDist, vF, vS, GlobalBlobID );
|
||||
ReorderBlobIDs(GlobalBlobID); // This will reorder by blob size
|
||||
printf("Identified %i blobs. Writing per-process output files. \n",nblobs);
|
||||
|
||||
double vF=0.0;
|
||||
double vS=0.0;
|
||||
double trimdist=1.0;
|
||||
printf("Execute blob identification algorithm... \n");
|
||||
// Loop over z=0 first -> blobs attached to this end considered "connected" for LB simulation
|
||||
i=0;
|
||||
int number=0;
|
||||
/* for (k=0;k<1;k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
if ( Phase(i,j,k) > vF ){
|
||||
if ( SignDist(i,j,k) > vS ){
|
||||
// node i,j,k is in the porespace
|
||||
number = number+ComputeBlob(blobs,nblobs,ncubes,GlobalBlobID,Phase,SignDist,vF,vS,i,j,k,temp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Specify the blob on the z axis
|
||||
if (ncubes > 0){
|
||||
b(nblobs) = number;
|
||||
// BlobList.push_back[number];
|
||||
printf("Number of non-wetting phase blobs is: %i \n",nblobs-1);
|
||||
nblobs++;
|
||||
}
|
||||
*/
|
||||
for (k=0;k<Nz;k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
for (i=0;i<Nx;i++){
|
||||
if ( GlobalBlobID(i,j,k) == -1 ){
|
||||
if ( Phase(i,j,k) > vF ){
|
||||
if ( SignDist(i,j,k) > vS ){
|
||||
// node i,j,k is in the porespace
|
||||
b(nblobs) = ComputeBlob(blobs,nblobs,ncubes,GlobalBlobID,Phase,SignDist,vF,vS,i,j,k,temp);
|
||||
nblobs++;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Otherwise, this point has already been assigned - ignore
|
||||
int sizeLoc = nx*ny*nz;
|
||||
int *LocalBlobID;
|
||||
LocalBlobID = new int [sizeLoc];
|
||||
|
||||
// Make sure list blob_nodes is large enough
|
||||
if ( nblobs > (int)b.length()-1){
|
||||
printf("Increasing size of blob list \n");
|
||||
b.resize(2*b.length());
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Go over all cubes again -> add any that do not contain nw phase
|
||||
bool add=1; // Set to false if any corners contain nw-phase ( F > vF)
|
||||
// int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners
|
||||
int count_in=0,count_out=0;
|
||||
int nodx,nody,nodz;
|
||||
for (k=0;k<Nz-1;k++){
|
||||
for (j=0;j<Ny-1;j++){
|
||||
for (i=0;i<Nx-1;i++){
|
||||
// Loop over cube corners
|
||||
add=1; // initialize to true - add unless corner occupied by nw-phase
|
||||
for (p=0;p<8;p++){
|
||||
nodx=i+cube[p][0];
|
||||
nody=j+cube[p][1];
|
||||
nodz=k+cube[p][2];
|
||||
if ( GlobalBlobID(nodx,nody,nodz) > -1 ){
|
||||
// corner occupied by nw-phase -> do not add
|
||||
add = 0;
|
||||
}
|
||||
}
|
||||
if ( add == 1 ){
|
||||
blobs(0,ncubes) = i;
|
||||
blobs(1,ncubes) = j;
|
||||
blobs(2,ncubes) = k;
|
||||
ncubes++;
|
||||
count_in++;
|
||||
}
|
||||
else { count_out++; }
|
||||
}
|
||||
}
|
||||
}
|
||||
b(nblobs) = count_in;
|
||||
nblobs++;
|
||||
printf("Identified %i blobs. Writing per-process output files. \n",nblobs);
|
||||
printf("File size (4 bytes per entry) %i, \n",sizeLoc);
|
||||
// read the files and populate main arrays
|
||||
for (int kproc=0; kproc<nprocz; kproc++){
|
||||
for (int jproc=0; jproc<nprocy; jproc++){
|
||||
for (int iproc=0; iproc<nprocx; iproc++){
|
||||
|
||||
int sizeLoc = nx*ny*nz;
|
||||
int *LocalBlobID;
|
||||
LocalBlobID = new int [sizeLoc];
|
||||
int proc = kproc*nprocx*nprocy + jproc*nprocx + iproc;
|
||||
|
||||
printf("File size (4 bytes per entry) %i, \n",sizeLoc);
|
||||
// read the files and populate main arrays
|
||||
for ( kproc=0; kproc<nprocz; kproc++){
|
||||
for ( jproc=0; jproc<nprocy; jproc++){
|
||||
for ( iproc=0; iproc<nprocx; iproc++){
|
||||
sprintf(LocalRankString,"%05d",proc);
|
||||
sprintf(LocalRankFilename,"%s%s","BlobLabel.",LocalRankString);
|
||||
|
||||
proc = kproc*nprocx*nprocy + jproc*nprocx + iproc;
|
||||
for (int k=0; k<nz; k++){
|
||||
for (int j=0; j<ny; j++){
|
||||
for (int i=0; i<nx; i++){
|
||||
//........................................................................
|
||||
int n = k*nx*ny+j*nx+i;
|
||||
//........................................................................
|
||||
int iglobal = iproc*(nx-2)+i;
|
||||
int jglobal = jproc*(ny-2)+j;
|
||||
int kglobal = kproc*(nz-2)+k;
|
||||
// periodic BC
|
||||
if (iglobal < 0 ) iglobal+=Nx;
|
||||
if (jglobal < 0 ) jglobal+=Ny;
|
||||
if (kglobal < 0 ) kglobal+=Nz;
|
||||
if (!(iglobal < Nx) ) iglobal-=Nx;
|
||||
if (!(jglobal < Ny) ) jglobal-=Ny;
|
||||
if (!(kglobal < Nz) ) kglobal-=Nz;
|
||||
//........................................................................
|
||||
LocalBlobID[n] = GlobalBlobID(iglobal,jglobal,kglobal);
|
||||
//........................................................................
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
sprintf(LocalRankString,"%05d",proc);
|
||||
sprintf(LocalRankFilename,"%s%s","BlobLabel.",LocalRankString);
|
||||
FILE *BLOBLOCAL = fopen(LocalRankFilename,"wb");
|
||||
fwrite(&LocalBlobID[0],4,sizeLoc,BLOBLOCAL);
|
||||
fclose(BLOBLOCAL);
|
||||
}
|
||||
}
|
||||
}
|
||||
printf("Wrote %i ranks of BlobLabel.xxxxx \n",nprocs);
|
||||
|
||||
for (k=0; k<nz; k++){
|
||||
for (j=0; j<ny; j++){
|
||||
for (i=0; i<nx; i++){
|
||||
//........................................................................
|
||||
n = k*nx*ny+j*nx+i;
|
||||
//........................................................................
|
||||
iglobal = iproc*(nx-2)+i;
|
||||
jglobal = jproc*(ny-2)+j;
|
||||
kglobal = kproc*(nz-2)+k;
|
||||
// periodic BC
|
||||
if (iglobal < 0 ) iglobal+=Nx;
|
||||
if (jglobal < 0 ) jglobal+=Ny;
|
||||
if (kglobal < 0 ) kglobal+=Nz;
|
||||
if (!(iglobal < Nx) ) iglobal-=Nx;
|
||||
if (!(jglobal < Ny) ) jglobal-=Ny;
|
||||
if (!(kglobal < Nz) ) kglobal-=Nz;
|
||||
//........................................................................
|
||||
LocalBlobID[n] = GlobalBlobID(iglobal,jglobal,kglobal);
|
||||
//........................................................................
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
FILE *BLOBLOCAL;
|
||||
BLOBLOCAL = fopen(LocalRankFilename,"wb");
|
||||
fwrite(&LocalBlobID[0],4,sizeLoc,BLOBLOCAL);
|
||||
fclose(BLOBLOCAL);
|
||||
}
|
||||
}
|
||||
}
|
||||
printf("Wrote %i ranks of BlobLabel.xxxxx \n",nprocs);
|
||||
|
||||
FILE *BLOBS;
|
||||
BLOBS = fopen("Blobs.dat","wb");
|
||||
fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
|
||||
fclose(BLOBS);
|
||||
|
||||
FILE *BLOBS = fopen("Blobs.dat","wb");
|
||||
fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
|
||||
fclose(BLOBS);
|
||||
MPI_Finalize();
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
|
131
tests/BlobIdentifyParallel.cpp
Normal file
131
tests/BlobIdentifyParallel.cpp
Normal file
|
@ -0,0 +1,131 @@
|
|||
// Sequential blob analysis
|
||||
// Reads parallel simulation data and performs connectivity analysis
|
||||
// and averaging on a blob-by-blob basis
|
||||
// James E. McClure 2014
|
||||
|
||||
#include <iostream>
|
||||
#include <math.h>
|
||||
#include "common/pmmc.h"
|
||||
#include "common/Communication.h"
|
||||
#include "analysis/analysis.h"
|
||||
#include "ProfilerApp.h"
|
||||
|
||||
|
||||
//#include "Domain.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
|
||||
inline void ReadBinaryFile(char *FILENAME, double *Data, int N)
|
||||
{
|
||||
int n;
|
||||
double value;
|
||||
ifstream File(FILENAME,ios::binary);
|
||||
for (n=0; n<N; n++){
|
||||
// Write the two density values
|
||||
File.read((char*) &value, sizeof(value));
|
||||
Data[n] = value;
|
||||
|
||||
}
|
||||
File.close();
|
||||
}
|
||||
|
||||
|
||||
void readRankData( int proc, int nx, int ny, int nz, DoubleArray& Phase, DoubleArray& SignDist )
|
||||
{
|
||||
Phase.resize(nx,ny,nz);
|
||||
SignDist.resize(nx,ny,nz);
|
||||
char file1[40], file2[40];
|
||||
sprintf(file1,"SignDist.%05d",proc);
|
||||
sprintf(file2,"Phase.%05d",proc);
|
||||
ReadBinaryFile(file1, Phase.get(), nx*ny*nz);
|
||||
ReadBinaryFile(file2, SignDist.get(), nx*ny*nz);
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
// Initialize MPI
|
||||
int rank, nprocs;
|
||||
MPI_Init(&argc,&argv);
|
||||
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
|
||||
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
|
||||
PROFILE_ENABLE(0);
|
||||
PROFILE_DISABLE_TRACE();
|
||||
PROFILE_SYNCHRONIZE();
|
||||
PROFILE_START("main");
|
||||
|
||||
if ( rank==0 ) {
|
||||
printf("-----------------------------------------------------------\n");
|
||||
printf("Labeling Blobs from Two-Phase Lattice Boltzmann Simulation \n");
|
||||
printf("-----------------------------------------------------------\n");
|
||||
}
|
||||
|
||||
//.......................................................................
|
||||
// Reading the domain information file
|
||||
//.......................................................................
|
||||
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
|
||||
double Lx, Ly, Lz;
|
||||
ifstream domain("Domain.in");
|
||||
domain >> nprocx;
|
||||
domain >> nprocy;
|
||||
domain >> nprocz;
|
||||
domain >> nx;
|
||||
domain >> ny;
|
||||
domain >> nz;
|
||||
domain >> nspheres;
|
||||
domain >> Lx;
|
||||
domain >> Ly;
|
||||
domain >> Lz;
|
||||
|
||||
// Check that the number of processors >= the number of ranks
|
||||
if ( rank==0 ) {
|
||||
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
|
||||
printf("Number of MPI ranks used: %i \n", nprocs);
|
||||
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
|
||||
}
|
||||
if ( nprocs < nprocx*nprocy*nprocz )
|
||||
ERROR("Insufficient number of processors");
|
||||
|
||||
// Get the rank info
|
||||
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
|
||||
|
||||
// Read the local file
|
||||
DoubleArray Phase;
|
||||
DoubleArray SignDist;
|
||||
readRankData( rank, nx+2, ny+2, nz+2, Phase, SignDist );
|
||||
|
||||
// Communication the halos
|
||||
fillHalo<double> fillData(rank_info,nx,ny,nz,1,1,1,0,1);
|
||||
fillData.fill(Phase);
|
||||
fillData.fill(SignDist);
|
||||
|
||||
// Find blob domains
|
||||
if ( rank==0 ) { printf("Finding blob domains\n"); }
|
||||
double vF=0.0;
|
||||
double vS=0.0;
|
||||
IntArray GlobalBlobID;
|
||||
int nblobs = ComputeGlobalBlobIDs(nx,ny,nz,rank_info,
|
||||
Phase,SignDist,vF,vS,GlobalBlobID);
|
||||
if ( rank==0 ) { printf("Identified %i blobs\n",nblobs); }
|
||||
|
||||
// Write the local blob ids
|
||||
char LocalRankFilename[100];
|
||||
sprintf(LocalRankFilename,"BlobLabel.%05i",rank);
|
||||
FILE *BLOBLOCAL = fopen(LocalRankFilename,"wb");
|
||||
fwrite(GlobalBlobID.get(),4,GlobalBlobID.length(),BLOBLOCAL);
|
||||
fclose(BLOBLOCAL);
|
||||
printf("Wrote BlobLabel.%05i \n",rank);
|
||||
|
||||
|
||||
/*FILE *BLOBS = fopen("Blobs.dat","wb");
|
||||
fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
|
||||
fclose(BLOBS);*/
|
||||
|
||||
PROFILE_STOP("main");
|
||||
PROFILE_SAVE("BlobIdentifyParallel",false);
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
MPI_Finalize();
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -8,6 +8,7 @@ INSTALL_LBPM_EXE( TestBubble )
|
|||
INSTALL_LBPM_EXE( BasicSimulator )
|
||||
INSTALL_LBPM_EXE( BlobAnalysis )
|
||||
INSTALL_LBPM_EXE( BlobIdentify )
|
||||
INSTALL_LBPM_EXE( BlobIdentifyParallel )
|
||||
|
||||
CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_DIR}/cylindertest COPYONLY )
|
||||
|
||||
|
|
|
@ -5,8 +5,9 @@
|
|||
#include <stdexcept>
|
||||
#include <fstream>
|
||||
|
||||
#include "Communication.h"
|
||||
#include "common/Communication.h"
|
||||
#include "common/MPI_Helpers.h"
|
||||
#include "common/Array.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
|
@ -182,11 +183,84 @@ int test_communication( MPI_Comm comm, int nprocx, int nprocy, int nprocz )
|
|||
}
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
int testHalo( MPI_Comm comm, int nprocx, int nprocy, int nprocz, int depth )
|
||||
{
|
||||
int rank,nprocs;
|
||||
MPI_Comm_rank(comm,&rank);
|
||||
MPI_Comm_size(comm,&nprocs);
|
||||
if ( rank==0 )
|
||||
printf("\nRunning Halo test %i %i %i %i\n",nprocx,nprocy,nprocz,depth);
|
||||
|
||||
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
|
||||
|
||||
int nx = 10;
|
||||
int ny = 11;
|
||||
int nz = 7;
|
||||
std::vector<size_t> size(4);
|
||||
size[0] = nx+2;
|
||||
size[1] = ny+2;
|
||||
size[2] = nz+2;
|
||||
size[3] = depth;
|
||||
Array<TYPE> array(size);
|
||||
array.fill(-1);
|
||||
|
||||
// Fill the local array
|
||||
int Nx = nx*nprocx;
|
||||
int Ny = ny*nprocy;
|
||||
int Nz = nz*nprocz;
|
||||
for (int i=0; i<nx; i++) {
|
||||
for (int j=0; j<ny; j++) {
|
||||
for (int k=0; k<nz; k++) {
|
||||
for (int d=0; d<depth; d++) {
|
||||
int iglobal = i + rank_info.ix*nx;
|
||||
int jglobal = j + rank_info.jy*ny;
|
||||
int kglobal = k + rank_info.kz*nz;
|
||||
int ijk = iglobal + jglobal*Nx + kglobal*Nx*Ny + d*Nx*Ny*Nz;
|
||||
array(i+1,j+1,k+1,d) = ijk;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Communicate the halo
|
||||
fillHalo<TYPE> fillData(rank_info,nx,ny,nz,1,1,1,0,depth);
|
||||
fillData.fill(array);
|
||||
|
||||
// Check the results
|
||||
bool pass = true;
|
||||
for (int i=-1; i<nx+1; i++) {
|
||||
for (int j=-1; j<ny+1; j++) {
|
||||
for (int k=-1; k<nz+1; k++) {
|
||||
for (int d=0; d<depth; d++) {
|
||||
int iglobal = i + rank_info.ix*nx;
|
||||
int jglobal = j + rank_info.jy*ny;
|
||||
int kglobal = k + rank_info.kz*nz;
|
||||
iglobal = (iglobal+Nx)%Nx;
|
||||
jglobal = (jglobal+Ny)%Ny;
|
||||
kglobal = (kglobal+Nz)%Nz;
|
||||
int ijk = iglobal + jglobal*Nx + kglobal*Nx*Ny + d*Nx*Ny*Nz;
|
||||
if ( array(i+1,j+1,k+1,d) != ijk )
|
||||
pass = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
int N_errors = 0;
|
||||
if ( !pass ) {
|
||||
std::cout << "Failed halo test\n";
|
||||
N_errors++;
|
||||
}
|
||||
return N_errors;
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
// Initialize MPI
|
||||
int nprocs;
|
||||
int rank,nprocs;
|
||||
MPI_Init(&argc,&argv);
|
||||
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
|
||||
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
|
||||
|
||||
// Run the test with different domains
|
||||
|
@ -200,11 +274,33 @@ int main(int argc, char **argv)
|
|||
N_errors += test_communication( MPI_COMM_WORLD, 1, 2, 2 );
|
||||
}
|
||||
|
||||
// Run the halo tests with different domains
|
||||
N_errors += testHalo<int>( MPI_COMM_WORLD, nprocs, 1, 1, 1 );
|
||||
N_errors += testHalo<int>( MPI_COMM_WORLD, 1, nprocs, 1, 1 );
|
||||
N_errors += testHalo<int>( MPI_COMM_WORLD, 1, 1, nprocs, 1 );
|
||||
N_errors += testHalo<double>( MPI_COMM_WORLD, nprocs, 1, 1, 3 );
|
||||
N_errors += testHalo<double>( MPI_COMM_WORLD, 1, nprocs, 1, 3 );
|
||||
N_errors += testHalo<double>( MPI_COMM_WORLD, 1, 1, nprocs, 3 );
|
||||
if ( nprocs==4 ) {
|
||||
N_errors += testHalo<int>( MPI_COMM_WORLD, 2, 2, 1, 1 );
|
||||
N_errors += testHalo<int>( MPI_COMM_WORLD, 2, 1, 2, 1 );
|
||||
N_errors += testHalo<int>( MPI_COMM_WORLD, 1, 2, 2, 1 );
|
||||
}
|
||||
if ( nprocs==8 ) {
|
||||
N_errors += testHalo<int>( MPI_COMM_WORLD, 2, 2, 2, 1 );
|
||||
}
|
||||
|
||||
// Finished
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
int N_errors_global=0;
|
||||
MPI_Allreduce( &N_errors, &N_errors_global, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
MPI_Finalize();
|
||||
if ( rank==0 ) {
|
||||
if ( N_errors_global==0 )
|
||||
std::cout << "All tests passed\n";
|
||||
else
|
||||
std::cout << "Some tests failed\n";
|
||||
}
|
||||
return N_errors_global;
|
||||
}
|
||||
|
|
Loading…
Reference in New Issue
Block a user