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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James McClure 2015-07-18 17:23:35 -04:00
commit c8b5e1c020
22 changed files with 1384 additions and 467 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;

2
IO/Writer.cpp
View File

@ -158,6 +158,7 @@ static std::vector<IO::MeshDatabase> writeMeshesNewFormat(
// Write the mesh data
void IO::writeData( int timestep, const std::vector<IO::MeshDataStruct>& meshData, int format )
{
PROFILE_START("writeData");
int rank = MPI_WORLD_RANK();
int size = MPI_WORLD_SIZE();
// Create the output directory
@ -196,6 +197,7 @@ void IO::writeData( int timestep, const std::vector<IO::MeshDataStruct>& meshDat
fprintf(fid,"%s\n",path);
fclose(fid);
}
PROFILE_STOP("writeData");
}

536
analysis/analysis.cpp
View File

@ -12,8 +12,8 @@ inline const TYPE* getPtr( const std::vector<TYPE>& x ) { return x.empty() ? NUL
/******************************************************************
* Compute the blobs *
******************************************************************/
int ComputePhaseComponent(IntArray &ComponentLabel,
IntArray &PhaseID, int VALUE, int &ncomponents,
static int ComputePhaseComponent( IntArray &ComponentLabel,
const IntArray &PhaseID, int VALUE, int &ncomponents,
int startx, int starty, int startz, IntArray &temp, bool periodic)
{
@ -309,7 +309,7 @@ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
return nblobs;
}
int ComputeLocalPhaseComponent(IntArray &PhaseID, int VALUE, IntArray &ComponentLabel, bool periodic )
int ComputeLocalPhaseComponent(const IntArray &PhaseID, int VALUE, IntArray &ComponentLabel, bool periodic )
{
PROFILE_START("ComputeLocalPhaseComponent");
size_t Nx = PhaseID.size(0);
@ -352,6 +352,7 @@ int ComputeLocalPhaseComponent(IntArray &PhaseID, int VALUE, IntArray &Component
return ncomponents;
}
/******************************************************************
* Reorder the global blob ids *
******************************************************************/
@ -412,6 +413,7 @@ void ReorderBlobIDs( IntArray& ID )
PROFILE_STOP("ReorderBlobIDs");
}
/******************************************************************
* Compute the global blob ids *
******************************************************************/
@ -468,7 +470,137 @@ static bool updateLocalIds( const std::map<int64_t,int64_t>& remote_map,
}
return changed;
}
int ComputeGlobalBlobIDs( int nx, int ny, int nz, RankInfoStruct rank_info,
static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info,
int nblobs, IntArray& IDs )
{
PROFILE_START("LocalToGlobalIDs",1);
const int rank = rank_info.rank[1][1][1];
int nprocs;
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
// Get the number of blobs for each rank
std::vector<int> N_blobs(nprocs,0);
PROFILE_START("LocalToGlobalIDs-Allgather",1);
MPI_Allgather(&nblobs,1,MPI_INT,getPtr(N_blobs),1,MPI_INT,MPI_COMM_WORLD);
PROFILE_STOP("LocalToGlobalIDs-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<IDs.length(); i++) {
if ( IDs(i) >= 0 )
IDs(i) += offset;
}
const Array<int> LocalIDs = IDs;
// Copy the ids and get the neighbors through the halos
fillHalo<int> fillData(rank_info,nx,ny,nz,1,1,1,0,1,true,true,true);
fillData.fill(IDs);
// 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)!=IDs(i) )
map[LocalIDs(i)].remote_ids.insert(IDs(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("LocalToGlobalIDs-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("LocalToGlobalIDs-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 = (IDs.size(0)-nx)/2;
int ngy = (IDs.size(1)-ny)/2;
int ngz = (IDs.size(2)-nz)/2;
for (size_t k=ngz; k<IDs.size(2)-ngz; k++) {
for (size_t j=ngy; j<IDs.size(1)-ngy; j++) {
for (size_t i=ngx; i<IDs.size(0)-ngx; i++) {
int id = IDs(i,j,k);
if ( id >= 0 )
IDs(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(IDs);
// Reorder based on size (and compress the id space
int N_blobs_global = ReorderBlobIDs2(IDs,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("LocalToGlobalIDs",1);
return N_blobs_global;
}
int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info,
const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS,
IntArray& GlobalBlobID )
{
@ -477,262 +609,176 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, RankInfoStruct rank_info,
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];
int nblobs = ComputeLocalBlobIDs(Phase,SignDist,vF,vS,GlobalBlobID,false);
// Compute the global ids
int nglobal = LocalToGlobalIDs( nx, ny, nz, rank_info, nblobs, GlobalBlobID );
PROFILE_STOP("ComputeGlobalBlobIDs");
return N_blobs_global;
return nglobal;
}
int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& rank_info,
const IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID )
{
PROFILE_START("ComputeGlobalPhaseComponent");
// First compute the local ids
int nblobs = ComputeLocalPhaseComponent(PhaseID,VALUE,GlobalBlobID,false);
// Compute the global ids
int nglobal = LocalToGlobalIDs( nx, ny, nz, rank_info, nblobs, GlobalBlobID );
PROFILE_STOP("ComputeGlobalPhaseComponent");
return nglobal;
}
int ComputeGlobalPhaseComponent( int nx, int ny, int nz, RankInfoStruct rank_info,
IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID )
/******************************************************************
* Compute the mapping of blob ids between timesteps *
******************************************************************/
void gatherSet( std::set<int>& set )
{
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 = ComputeLocalPhaseComponent(PhaseID,VALUE,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;
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);
}
// 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));
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));
}
}
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));
// 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);
}
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;
// 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);
}
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];
// 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));
}
}
}
// 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;
// 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;
}

41
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,
int ComputeLocalPhaseComponent( const 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
@ -68,14 +73,15 @@ int ComputeLocalPhaseComponent(IntArray &PhaseID, int VALUE, IntArray &Component
* @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,
int ComputeGlobalBlobIDs( int nx, int ny, int nz, const 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
* @return Returns the number of cubes in the 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
@ -85,9 +91,9 @@ 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, const RankInfoStruct& rank_info,
const IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID );
int ComputeGlobalPhaseComponent( int nx, int ny, int nz, RankInfoStruct rank_info,
IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID );
/*!
* @brief Reorder the blobs
@ -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

12
common/Communication.hpp
View File

@ -133,11 +133,11 @@ void fillHalo<TYPE>::fill( Array<TYPE>& data )
}
}
}
// Recv the dst data and unpack
// Recv the dst data and unpack (we recive in reverse order to match the sends)
MPI_Status status;
for (int i=0; i<3; i++) {
for (int j=0; j<3; j++) {
for (int k=0; k<3; k++) {
for (int i=2; i>=0; i--) {
for (int j=2; j>=0; j--) {
for (int k=2; k>=0; k--) {
if ( !fill_pattern[i][j][k] )
continue;
MPI_Wait(&recv_req[i][j][k],&status);
@ -206,6 +206,7 @@ template<class TYPE>
template<class TYPE1, class TYPE2>
void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
{
PROFILE_START("fillHalo::copy",1);
ASSERT(src.size(0)==nx||src.size(0)==nx+2*ngx);
ASSERT(dst.size(0)==nx||dst.size(0)==nx+2*ngx);
bool src_halo = src.size(0)==nx+2*ngx;
@ -242,7 +243,7 @@ void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
}
}
} else if ( !src_halo && dst_halo ) {
// Src has halos
// Dst has halos
for (size_t k=0; k<nz; k++) {
for (size_t j=0; j<ny; j++) {
for (size_t i=0; i<nx; i++) {
@ -252,6 +253,7 @@ void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
}
fill(dst);
}
PROFILE_STOP("fillHalo::copy",1);
}

8
common/MPI_Helpers.cpp
View File

@ -168,7 +168,6 @@ int MPI_Wait( MPI_Request*, MPI_Status* )
}
int MPI_Bcast( void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm )
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Send(const void *buf, int count, MPI_Datatype datatype, int dest, int tag,
@ -208,6 +207,13 @@ int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
ERROR("Not implimented yet");
return 0;
}
int MPI_Allgatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
void *recvbuf, const int *recvcounts, const int *displs,
MPI_Datatype recvtype, MPI_Comm comm)
{
ERROR("Not implimented yet");
return 0;
}
int MPI_Sendrecv(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
int dest, int sendtag,
void *recvbuf, int recvcount, MPI_Datatype recvtype,

3
common/MPI_Helpers.h
View File

@ -44,6 +44,9 @@
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
void *recvbuf, int recvcount, MPI_Datatype recvtype,
MPI_Comm comm);
int MPI_Allgatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
void *recvbuf, const int *recvcounts, const int *displs,
MPI_Datatype recvtype, MPI_Comm comm);
int MPI_Sendrecv(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
int dest, int sendtag,
void *recvbuf, int recvcount, MPI_Datatype recvtype,

6
common/ScaLBL.h
View File

@ -85,6 +85,12 @@ extern "C" void ColorBC_inlet(double *Phi, double *Den, double *A_even, double *
extern "C" void ColorBC_outlet(double *Phi, double *Den, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz);
extern "C" void ScaLBL_D3Q19_Velocity_BC_z(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz);
extern "C" void ScaLBL_D3Q19_Velocity_BC_Z(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz, int outlet);
class ScaLBL_Communicator{
public:
//......................................................................................

388
common/StackTrace.cpp Normal file
View File

@ -0,0 +1,388 @@
#include "StackTrace.h"
#include <iostream>
#include <sstream>
#include <cstring>
#include <algorithm>
#if __cplusplus > 199711L
#include <mutex>
#endif
// Detect the OS and include system dependent headers
#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64) || defined(_MSC_VER)
// Note: windows has not been testeds
#define USE_WINDOWS
#include <iostream>
#include <windows.h>
#include <process.h>
#include <stdio.h>
#include <tchar.h>
#include <psapi.h>
#include <DbgHelp.h>
//#pragma comment(lib, psapi.lib) //added
//#pragma comment(linker, /DEFAULTLIB:psapi.lib)
#elif defined(__APPLE__)
#define USE_MAC
#include <sys/time.h>
#include <sys/sysctl.h>
#include <signal.h>
#include <execinfo.h>
#include <dlfcn.h>
#include <mach/mach.h>
#include <sys/types.h>
#include <sys/sysctl.h>
#include <unistd.h>
#include <sched.h>
#elif defined(__linux) || defined(__unix) || defined(__posix)
#define USE_LINUX
#define USE_NM
#include <sys/time.h>
#include <time.h>
#include <execinfo.h>
#include <dlfcn.h>
#include <malloc.h>
#include <unistd.h>
#include <sched.h>
#else
#error Unknown OS
#endif
#ifdef __GNUC__
#define USE_ABI
#include <cxxabi.h>
#endif
#ifndef NULL_USE
#define NULL_USE(variable) do { \
if(0) {char *temp = (char *)&variable; temp++;} \
}while(0)
#endif
// Utility to strip the path from a filename
inline std::string stripPath( const std::string& filename )
{
if ( filename.empty() ) { return std::string(); }
int i=0;
for (i=(int)filename.size()-1; i>=0&&filename[i]!=47&&filename[i]!=92; i--) {}
i = std::max(0,i+1);
return filename.substr(i);
}
// Inline function to subtract two addresses returning the absolute difference
inline void* subtractAddress( void* a, void* b ) {
return reinterpret_cast<void*>( std::abs(
reinterpret_cast<long long int>(a)-reinterpret_cast<long long int>(b) ) );
}
/****************************************************************************
* stack_info *
****************************************************************************/
std::string StackTrace::stack_info::print() const
{
char tmp[32];
sprintf(tmp,"0x%016llx: ",reinterpret_cast<unsigned long long int>(address));
std::string stack(tmp);
sprintf(tmp,"%i",line);
std::string line_str(tmp);
stack += stripPath(object);
stack.resize(std::max<size_t>(stack.size(),38),' ');
stack += " " + function;
if ( !filename.empty() && line>0 ) {
stack.resize(std::max<size_t>(stack.size(),70),' ');
stack += " " + stripPath(filename) + ":" + line_str;
} else if ( !filename.empty() ) {
stack.resize(std::max<size_t>(stack.size(),70),' ');
stack += " " + stripPath(filename);
} else if ( line>0 ) {
stack += " : " + line_str;
}
return stack;
}
/****************************************************************************
* Function to find an entry *
****************************************************************************/
template <class TYPE>
inline size_t findfirst( const std::vector<TYPE>& X, TYPE Y )
{
if ( X.empty() )
return 0;
size_t lower = 0;
size_t upper = X.size()-1;
if ( X[lower] >= Y )
return lower;
if ( X[upper] < Y )
return upper;
while ( (upper-lower) != 1 ) {
size_t value = (upper+lower)/2;
if ( X[value] >= Y )
upper = value;
else
lower = value;
}
return upper;
}
/****************************************************************************
* Function to get symbols for the executable from nm (if availible) *
* Note: this function maintains an internal cached copy to prevent *
* exccessive calls to nm. This function also uses a lock to ensure *
* thread safety. *
****************************************************************************/
#if __cplusplus <= 199711L
class mutex_class {
public:
void lock() {}
void unlock() {}
};
mutex_class getSymbols_mutex;
#else
std::mutex getSymbols_mutex;
#endif
struct global_symbols_struct {
std::vector<void*> address;
std::vector<char> type;
std::vector<std::string> obj;
int error;
} global_symbols;
static std::string get_executable()
{
std::string exe;
try {
#ifdef USE_LINUX
char *buf = new char[0x10000];
int len = ::readlink("/proc/self/exe",buf,0x10000);
if ( len!=-1 ) {
buf[len] = '\0';
exe = std::string(buf);
}
#endif
} catch (...) {}
return exe;
}
std::string global_exe_name = get_executable();
static const global_symbols_struct& getSymbols2( )
{
static bool loaded = false;
static global_symbols_struct data;
// Load the symbol tables if they have not been loaded
if ( !loaded ) {
getSymbols_mutex.lock();
if ( !loaded ) {
loaded = true;
#ifdef USE_NM
try {
char cmd[1024];
sprintf(cmd,"nm --demangle --numeric-sort %s",global_exe_name.c_str());
FILE *in = popen(cmd,"r");
if ( in==NULL ) {
data.error = -2;
return data;
}
char *buf = new char[0x100000];
while ( fgets(buf,0xFFFFF,in)!=NULL ) {
if ( buf[0]==' ' || buf==NULL )
continue;
char *a = buf;
char *b = strchr(a,' '); if (b==NULL) {continue;} b[0] = 0; b++;
char *c = strchr(b,' '); if (c==NULL) {continue;} c[0] = 0; c++;
char *d = strchr(c,'\n'); if ( d ) { d[0]=0; }
size_t add = strtoul(a,NULL,16);
data.address.push_back( reinterpret_cast<void*>(add) );
data.type.push_back( b[0] );
data.obj.push_back( std::string(c) );
}
pclose(in);
delete [] buf;
} catch (...) {
data.error = -3;
}
data.error = 0;
#else
data.error = -1;
#endif
}
getSymbols_mutex.unlock();
}
return data;
}
int StackTrace::getSymbols( std::vector<void*>& address, std::vector<char>& type,
std::vector<std::string>& obj )
{
const global_symbols_struct& data = getSymbols2();
address = data.address;
type = data.type;
obj = data.obj;
return data.error;
}
/****************************************************************************
* Function to get the current call stack *
****************************************************************************/
static void getFileAndLine( StackTrace::stack_info& info )
{
#if defined(USE_LINUX) || defined(USE_MAC)
void *address = info.address;
if ( info.object.find(".so")!=std::string::npos )
address = info.address2;
char buf[4096];
sprintf(buf, "addr2line -C -e %s -f -i %lx 2> /dev/null",
info.object.c_str(),reinterpret_cast<unsigned long int>(address));
FILE* f = popen(buf, "r");
if (f == NULL)
return;
buf[4095] = 0;
// get function name
char *rtn = fgets(buf,4095,f);
if ( info.function.empty() && rtn==buf ) {
info.function = std::string(buf);
info.function.resize(std::max<size_t>(info.function.size(),1)-1);
}
// get file and line
rtn = fgets(buf,4095,f);
if ( buf[0]!='?' && buf[0]!=0 && rtn==buf ) {
size_t i = 0;
for (i=0; i<4095 && buf[i]!=':'; i++) { }
info.filename = std::string(buf,i);
info.line = atoi(&buf[i+1]);
}
pclose(f);
#endif
}
// Try to use the global symbols to decode info about the stack
static void getDataFromGlobalSymbols( StackTrace::stack_info& info )
{
const global_symbols_struct& data = getSymbols2();
if ( data.error==0 ) {
size_t index = findfirst(global_symbols.address,info.address);
if ( index > 0 )
info.object = global_symbols.obj[index-1];
else
info.object = global_exe_name;
}
}
StackTrace::stack_info StackTrace::getStackInfo( void* address )
{
StackTrace::stack_info info;
info.address = address;
#ifdef _GNU_SOURCE
Dl_info dlinfo;
if ( !dladdr(address, &dlinfo) ) {
getDataFromGlobalSymbols( info );
getFileAndLine(info);
return info;
}
info.address2 = subtractAddress(info.address,dlinfo.dli_fbase);
info.object = std::string(dlinfo.dli_fname);
#if defined(USE_ABI)
int status;
char *demangled = abi::__cxa_demangle(dlinfo.dli_sname,NULL,0,&status);
if ( status == 0 && demangled!=NULL ) {
info.function = std::string(demangled);
} else if ( dlinfo.dli_sname!=NULL ) {
info.function = std::string(dlinfo.dli_sname);
}
free(demangled);
#else
if ( dlinfo.dli_sname!=NULL )
info.function = std::string(dlinfo.dli_sname);
#endif
#else
getDataFromGlobalSymbols( info );
#endif
// Get the filename / line number
getFileAndLine(info);
return info;
}
std::vector<StackTrace::stack_info> StackTrace::getCallStack()
{
std::vector<StackTrace::stack_info> stack_list;
#if defined(USE_LINUX) || defined(USE_MAC)
// Get the trace
void *trace[100];
memset(trace,0,100*sizeof(void*));
int trace_size = backtrace(trace,100);
stack_list.reserve(trace_size);
for (int i=0; i<trace_size; ++i)
stack_list.push_back(getStackInfo(trace[i]));
#elif defined(USE_WINDOWS)
#ifdef DBGHELP
::CONTEXT lContext;
::ZeroMemory( &lContext, sizeof( ::CONTEXT ) );
::RtlCaptureContext( &lContext );
::STACKFRAME64 lFrameStack;
::ZeroMemory( &lFrameStack, sizeof( ::STACKFRAME64 ) );
lFrameStack.AddrPC.Offset = lContext.Rip;
lFrameStack.AddrFrame.Offset = lContext.Rbp;
lFrameStack.AddrStack.Offset = lContext.Rsp;
lFrameStack.AddrPC.Mode = lFrameStack.AddrFrame.Mode = lFrameStack.AddrStack.Mode = AddrModeFlat;
#ifdef _M_IX86
DWORD MachineType = IMAGE_FILE_MACHINE_I386;
#endif
#ifdef _M_X64
DWORD MachineType = IMAGE_FILE_MACHINE_AMD64;
#endif
#ifdef _M_IA64
DWORD MachineType = IMAGE_FILE_MACHINE_IA64;
#endif
while ( 1 ) {
int rtn = ::StackWalk64( MachineType, ::GetCurrentProcess(), ::GetCurrentThread(),
&lFrameStack, MachineType == IMAGE_FILE_MACHINE_I386 ? 0 : &lContext,
NULL, &::SymFunctionTableAccess64, &::SymGetModuleBase64, NULL );
if( !rtn )
break;
if( lFrameStack.AddrPC.Offset == 0 )
break;
::MEMORY_BASIC_INFORMATION lInfoMemory;
::VirtualQuery( ( ::PVOID )lFrameStack.AddrPC.Offset, &lInfoMemory, sizeof( lInfoMemory ) );
if ( lInfoMemory.Type==MEM_PRIVATE )
continue;
::DWORD64 lBaseAllocation = reinterpret_cast< ::DWORD64 >( lInfoMemory.AllocationBase );
::TCHAR lNameModule[ 1024 ];
::HMODULE hBaseAllocation = reinterpret_cast< ::HMODULE >( lBaseAllocation );
::GetModuleFileName( hBaseAllocation, lNameModule, 1024 );
PIMAGE_DOS_HEADER lHeaderDOS = reinterpret_cast<PIMAGE_DOS_HEADER>( lBaseAllocation );
if ( lHeaderDOS==NULL )
continue;
PIMAGE_NT_HEADERS lHeaderNT = reinterpret_cast<PIMAGE_NT_HEADERS>( lBaseAllocation + lHeaderDOS->e_lfanew );
PIMAGE_SECTION_HEADER lHeaderSection = IMAGE_FIRST_SECTION( lHeaderNT );
::DWORD64 lRVA = lFrameStack.AddrPC.Offset - lBaseAllocation;
::DWORD64 lNumberSection = ::DWORD64();
::DWORD64 lOffsetSection = ::DWORD64();
for( int lCnt = ::DWORD64(); lCnt < lHeaderNT->FileHeader.NumberOfSections; lCnt++, lHeaderSection++ ) {
::DWORD64 lSectionBase = lHeaderSection->VirtualAddress;
::DWORD64 lSectionEnd = lSectionBase + std::max<::DWORD64>(
lHeaderSection->SizeOfRawData, lHeaderSection->Misc.VirtualSize );
if( ( lRVA >= lSectionBase ) && ( lRVA <= lSectionEnd ) ) {
lNumberSection = lCnt + 1;
lOffsetSection = lRVA - lSectionBase;
//break;
}
}
StackTrace::stack_info info;
info.object = lNameModule;
info.address = reinterpret_cast<void*>(lRVA);
char tmp[20];
sprintf(tmp,"0x%016llx",static_cast<unsigned long long int>(lOffsetSection));
info.function = std::to_string(lNumberSection) + ":" + std::string(tmp);
stack_list.push_back(info);
}
#endif
#else
#warning Stack trace is not supported on this compiler/OS
#endif
return stack_list;
}

47
common/StackTrace.h Normal file
View File

@ -0,0 +1,47 @@
#ifndef included_StackTrace
#define included_StackTrace
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <vector>
namespace StackTrace {
struct stack_info {
void *address;
void *address2;
std::string object;
std::string function;
std::string filename;
int line;
//! Default constructor
stack_info(): address(NULL), address2(NULL), line(0) {}
//! Print the stack info
std::string print() const;
};
//! Function to return the current call stack
std::vector<stack_info> getCallStack();
//! Function to return the stack info for a given address
stack_info getStackInfo( void* address );
/*!
* Return the symbols from the current executable (not availible for all platforms)
* @return Returns 0 if sucessful
*/
int getSymbols( std::vector<void*>& address, std::vector<char>& type, std::vector<std::string>& obj );
} // namespace StackTrace
#endif

115
common/TwoPhase.h
View File

@ -424,7 +424,8 @@ void TwoPhase::UpdateMeshValues(){
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
if (SDs(i,j,k) < 0.0){
n = k*Nx*Ny+j*Nx+i;
if (Dm.id[n] == 0){
// Solid phase
PhaseID(i,j,k) = 0;
}
@ -1007,64 +1008,68 @@ void TwoPhase::PrintComponents(int timestep){
if (Dm.rank==0){
printf("PRINT COMPONENT AVEREAGES: time = %i \n",timestep);
for (int b=0; b<NumberComponents_NWP; b++){
fprintf(NWPLOG,"%i ",timestep-5);
fprintf(NWPLOG,"%i ",b);
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VOL,b));
// fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(TRIMVOL,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(PRS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(AWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(ANS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(JWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(KWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(LWNS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CWNS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNSX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNSY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNSZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VSQ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNXX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNYY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNZZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNXY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNXZ,b));
fprintf(NWPLOG,"%.5g\n",ComponentAverages_NWP(GWNYZ,b));
if (ComponentAverages_NWP(TRIMVOL,b) > 0.0){
fprintf(NWPLOG,"%i ",timestep-5);
fprintf(NWPLOG,"%i ",b);
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VOL,b));
// fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(TRIMVOL,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(PRS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(AWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(ANS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(JWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(KWN,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(LWNS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CWNS,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNSX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNSY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VWNSZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(VSQ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNXX,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNYY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNZZ,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNXY,b));
fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(GWNXZ,b));
fprintf(NWPLOG,"%.5g\n",ComponentAverages_NWP(GWNYZ,b));
}
}
fflush(NWPLOG);
for (int b=0; b<NumberComponents_WP; b++){
fprintf(WPLOG,"%i ",timestep-5);
fprintf(WPLOG,"%i ",b);
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VOL,b));
// fprintf(WPLOG,"%.5g ",ComponentAverages_WP(TRIMVOL,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(PRS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(AWN,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(AWS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(JWN,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(KWN,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(LWNS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CWNS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNSX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNSY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNSZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VSQ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNXX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNYY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNZZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNXY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNXZ,b));
fprintf(WPLOG,"%.5g\n",ComponentAverages_WP(GWNYZ,b));
if (ComponentAverages_WP(TRIMVOL,b) > 0.0){
fprintf(WPLOG,"%i ",timestep-5);
fprintf(WPLOG,"%i ",b);
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VOL,b));
// fprintf(WPLOG,"%.5g ",ComponentAverages_WP(TRIMVOL,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(PRS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(AWN,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(AWS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(JWN,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(KWN,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(LWNS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CWNS,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNSX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNSY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VWNSZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(VSQ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNXX,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNYY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNZZ,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNXY,b));
fprintf(WPLOG,"%.5g ",ComponentAverages_WP(GWNXZ,b));
fprintf(WPLOG,"%.5g\n",ComponentAverages_WP(GWNYZ,b));
}
}
fflush(WPLOG);

120
common/Utilities.cpp
View File

@ -1,5 +1,5 @@
#include "common/Utilities.h"
#include "common/StackTrace.h"
#include <iostream>
#include <sstream>
@ -76,11 +76,11 @@ void Utilities::abort(const std::string &message, const std::string &filename, c
msg << "Bytes used = " << N_bytes << std::endl;
}
if ( abort_printStack ) {
std::vector<std::string> stack = Utilities::getCallStack();
std::vector<StackTrace::stack_info> stack = StackTrace::getCallStack();
msg << std::endl;
msg << "Stack Trace:\n";
for (size_t i=0; i<stack.size(); i++)
msg << " " << stack[i] << std::endl;
msg << " " << stack[i].print() << std::endl;
}
msg << std::endl << message << std::endl;
// Print the message and abort
@ -132,28 +132,33 @@ void MPI_error_handler_fun( MPI_Comm *comm, int *err, ... )
/****************************************************************************
* Function to handle unhandled exceptions *
****************************************************************************/
bool tried_MPI_Abort=false;
void term_func_abort(int err)
{
printf("Exiting due to abort (%i)\n",err);
std::vector<std::string> stack = Utilities::getCallStack();
std::vector<StackTrace::stack_info> stack = StackTrace::getCallStack();
std::string message = "Stack Trace:\n";
for (size_t i=0; i<stack.size(); i++)
message += " " + stack[i] += "\n";
message += " " + stack[i].print() += "\n";
message += "\nExiting\n";
// Print the message and abort
std::cerr << message;
#ifdef USE_MPI
if ( !abort_throwException )
if ( !abort_throwException && !tried_MPI_Abort ) {
tried_MPI_Abort = true;
MPI_Abort(MPI_COMM_WORLD,-1);
}
#endif
exit(-1);
}
static int tried_throw = 0;
#if defined(USE_LINUX) || defined(USE_MAC)
static int tried_throw = 0;
#endif
void term_func()
{
// Try to re-throw the last error to get the last message
std::string last_message;
#ifdef USE_LINUX
#if defined(USE_LINUX) || defined(USE_MAC)
try {
if ( tried_throw==0 ) {
tried_throw = 1;
@ -253,105 +258,8 @@ size_t Utilities::getMemoryUsage()
/****************************************************************************
* Function to get the current call stack *
* Functions to get the time and timer resolution *
****************************************************************************/
std::vector<std::string> Utilities::getCallStack()
{
std::vector<std::string> stack_list;
#if defined(USE_ABI)
void *trace[100];
memset(trace,0,100*sizeof(void*));
Dl_info dlinfo;
int status;
const char *symname;
char *demangled=NULL;
int trace_size = backtrace(trace,100);
for (int i=0; i<trace_size; ++i) {
if(!dladdr(trace[i], &dlinfo))
continue;
symname = dlinfo.dli_sname;
demangled = abi::__cxa_demangle(symname, NULL, 0, &status);
if(status == 0 && demangled)
symname = demangled;
std::string object = std::string(dlinfo.dli_fname);
std::string function = "";
if ( symname!=NULL )
function = std::string(symname);
if ( i!=0 ) { // Skip the current function
std::string stack_item = object + ": " + function;
//stack_item = "object: " + object;
//stack_item += "function: " + function;
stack_list.push_back(stack_item);
}
if ( demangled!=NULL ) {
free(demangled);
demangled=NULL;
}
}
#elif defined(USE_WINDOWS)
::CONTEXT lContext;
::ZeroMemory( &lContext, sizeof( ::CONTEXT ) );
::RtlCaptureContext( &lContext );
::STACKFRAME64 lFrameStack;
::ZeroMemory( &lFrameStack, sizeof( ::STACKFRAME64 ) );
lFrameStack.AddrPC.Offset = lContext.Rip;
lFrameStack.AddrFrame.Offset = lContext.Rbp;
lFrameStack.AddrStack.Offset = lContext.Rsp;
lFrameStack.AddrPC.Mode = lFrameStack.AddrFrame.Mode = lFrameStack.AddrStack.Mode = AddrModeFlat;
#ifdef _M_IX86
DWORD MachineType = IMAGE_FILE_MACHINE_I386;
#endif
#ifdef _M_X64
DWORD MachineType = IMAGE_FILE_MACHINE_AMD64;
#endif
#ifdef _M_IA64
DWORD MachineType = IMAGE_FILE_MACHINE_IA64;
#endif
while ( 1 ) {
int rtn = ::StackWalk64( MachineType, ::GetCurrentProcess(), ::GetCurrentThread(),
&lFrameStack, MachineType == IMAGE_FILE_MACHINE_I386 ? 0 : &lContext,
NULL, &::SymFunctionTableAccess64, &::SymGetModuleBase64, NULL );
if( !rtn )
break;
if( lFrameStack.AddrPC.Offset == 0 )
break;
::MEMORY_BASIC_INFORMATION lInfoMemory;
::VirtualQuery( ( ::PVOID )lFrameStack.AddrPC.Offset, &lInfoMemory, sizeof( lInfoMemory ) );
if ( lInfoMemory.Type==MEM_PRIVATE )
continue;
::DWORD64 lBaseAllocation = reinterpret_cast< ::DWORD64 >( lInfoMemory.AllocationBase );
::TCHAR lNameModule[ 1024 ];
::HMODULE hBaseAllocation = reinterpret_cast< ::HMODULE >( lBaseAllocation );
::GetModuleFileName( hBaseAllocation, lNameModule, 1024 );
PIMAGE_DOS_HEADER lHeaderDOS = reinterpret_cast<PIMAGE_DOS_HEADER>( lBaseAllocation );
if ( lHeaderDOS==NULL )
continue;
PIMAGE_NT_HEADERS lHeaderNT = reinterpret_cast<PIMAGE_NT_HEADERS>( lBaseAllocation + lHeaderDOS->e_lfanew );
PIMAGE_SECTION_HEADER lHeaderSection = IMAGE_FIRST_SECTION( lHeaderNT );
::DWORD64 lRVA = lFrameStack.AddrPC.Offset - lBaseAllocation;
::DWORD64 lNumberSection = ::DWORD64();
::DWORD64 lOffsetSection = ::DWORD64();
for( int lCnt = ::DWORD64(); lCnt < lHeaderNT->FileHeader.NumberOfSections; lCnt++, lHeaderSection++ ) {
::DWORD64 lSectionBase = lHeaderSection->VirtualAddress;
::DWORD64 lSectionEnd = lSectionBase + max( lHeaderSection->SizeOfRawData, lHeaderSection->Misc.VirtualSize );
if( ( lRVA >= lSectionBase ) && ( lRVA <= lSectionEnd ) ) {
lNumberSection = lCnt + 1;
lOffsetSection = lRVA - lSectionBase;
break;
}
}
std::stringstream stream;
stream << lNameModule << " : 000" << lNumberSection << " : " << reinterpret_cast<void*>(lOffsetSection);
stack_list.push_back(stream.str());
}
#else
#warning Stack trace is not supported on this compiler/OS
#endif
return stack_list;
}
// Functions to get the time and timer resolution
#if defined(USE_WINDOWS)
double Utilities::time()
{

2
common/Utilities.h
View File

@ -51,8 +51,6 @@ namespace Utilities
*/
size_t getMemoryUsage();
//! Function to get the current call stack
std::vector<std::string> getCallStack();
//! Function to get an arbitrary point in time
double time();

1
cpu/Color.cpp
View File

@ -336,6 +336,7 @@ extern "C" void PressureBC_outlet(double *disteven, double *distodd, double dout
}
}
//*************************************************************************
extern "C" void ComputeColorGradient(char *ID, double *phi, double *ColorGrad, int Nx, int Ny, int Nz)
{

113
cpu/D3Q19.cpp
View File

@ -227,6 +227,119 @@ extern "C" void SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, i
}
}
extern "C" void ScaLBL_D3Q19_Velocity_BC_z(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz)
{
int n,N;
// distributions
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double din;
N = Nx*Ny*Nz;
for (n=Nx*Ny; n<2*Nx*Ny; n++){
//........................................................................
// Read distributions from "opposite" memory convention
//........................................................................
//........................................................................
f1 = distodd[n];
f3 = distodd[N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//........................................................................
f0 = disteven[n];
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//...................................................
// Determine the outlet flow velocity
// uz = 1.0 - (f0+f4+f3+f2+f1+f8+f7+f9+f10 +
// 2*(f5+f15+f18+f11+f14))/din;
din = (f0+f4+f3+f2+f1+f8+f7+f9+f10+2*(f5+f15+f18+f11+f14))/(1.0-uz);
// Set the unknown distributions:
f6 = f5 + 0.3333333333333333*din*uz;
f16 = f15 + 0.1666666666666667*din*uz;
f17 = f16 + f4 - f3-f15+f18+f8-f7 +f9-f10;
f12= (din*uz+f5+ f15+f18+f11+f14-f6-f16-f17-f2+f1-f14+f11-f8+f7+f9-f10)*0.5;
f13= din*uz+f5+ f15+f18+f11+f14-f6-f16-f17-f12;
//........Store in "opposite" memory location..........
disteven[3*N+n] = f6;
disteven[6*N+n] = f12;
distodd[6*N+n] = f13;
disteven[8*N+n] = f16;
distodd[8*N+n] = f17;
//...................................................
}
}
extern "C" void ScaLBL_D3Q19_Velocity_BC_Z(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz, int outlet)
{
int n,N;
// distributions
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double dout;
N = Nx*Ny*Nz;
// Loop over the boundary - threadblocks delineated by start...finish
for (n=outlet; n<N-Nx*Ny; n++){
//........................................................................
// Read distributions from "opposite" memory convention
//........................................................................
f1 = distodd[n];
f3 = distodd[N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//........................................................................
f0 = disteven[n];
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//uz = -1.0 + (f0+f4+f3+f2+f1+f8+f7+f9+f10 + 2*(f6+f16+f17+f12+f13))/dout;
dout = (f0+f4+f3+f2+f1+f8+f7+f9+f10 + 2*(f6+f16+f17+f12+f13))/(1.0+uz);
f5 = f6 - 0.33333333333333338*dout* uz;
f15 = f16 - 0.16666666666666678*dout* uz;
f18 = f15 - f4 + f3-f16+f17-f8+f7-f9+f10;
f11 = (-dout*uz+f6+ f16+f17+f12+f13-f5-f15-f18+f2-f1-f13+f12+f8-f7-f9+f10)*0.5;
f14 = -dout*uz+f6+ f16+f17+f12+f13-f5-f15-f18-f11;
//........Store in "opposite" memory location..........
distodd[2*N+n] = f5;
distodd[5*N+n] = f11;
disteven[7*N+n] = f14;
distodd[7*N+n] = f15;
disteven[9*N+n] = f18;
//...................................................
}
}
extern "C" void ComputeVelocityD3Q19(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz)
{
int n,N;

25
example/MicroModel/Color.in.pressures Normal file
View File

@ -0,0 +1,25 @@
0 1 1.0053563 0.9946437
0 1 1.0078242 0.9921758
0 1 1.0091379 0.9908621
0 1 1.0109272 0.9890728
0 1 1.0128383 0.9871617
0 1 1.0153004 0.9846996
0 1 1.0165677 0.9834323
0 1 1.0176581 0.9823419
0 1 1.0181366 0.9818634
0 1 1.0184295 0.9815705
0 1 1.0187746 0.9812254
0 1 1.0190356 0.9809644
0 1 1.0193633 0.9806367
0 1 1.0195025 0.9804975
0 1 1.0200448 0.9799552
0 1 1.0202942 0.9797058
0 1 1.0215934 0.9784066
0 1 1.0227563 0.9772437
0 1 1.0246442 0.9753558
0 1 1.0264103 0.9735897
0 1 1.0289913 0.9710087
0 1 1.0315926 0.9684074
0 1 1.0339619 0.9660381
0 1 1.0401737 0.9598263
0 1 1.046342 0.953658

124
gpu/D3Q19.cu
View File

@ -349,7 +349,117 @@ __global__ void dvc_ComputePressureD3Q19(char *ID, double *disteven, double *di
}
}
__global__ void dvc_D3Q19_Velocity_BC_z(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz)
{
int n,N;
// distributions
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double uz;
N = Nx*Ny*Nz;
n = Nx*Ny + blockIdx.x*blockDim.x + threadIdx.x;
if (n < 2*Nx*Ny){
//........................................................................
// Read distributions from "opposite" memory convention
//........................................................................
//........................................................................
f1 = distodd[n];
f3 = distodd[N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//........................................................................
f0 = disteven[n];
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//...................................................
// Determine the outlet flow velocity
// uz = 1.0 - (f0+f4+f3+f2+f1+f8+f7+f9+f10 +
// 2*(f5+f15+f18+f11+f14))/din;
din = (f0+f4+f3+f2+f1+f8+f7+f9+f10+2*(f5+f15+f18+f11+f14))/(1.0-uz);
// Set the unknown distributions:
f6 = f5 + 0.3333333333333333*din*uz;
f16 = f15 + 0.1666666666666667*din*uz;
f17 = f16 + f4 - f3-f15+f18+f8-f7 +f9-f10;
f12= (din*uz+f5+ f15+f18+f11+f14-f6-f16-f17-f2+f1-f14+f11-f8+f7+f9-f10)*0.5;
f13= din*uz+f5+ f15+f18+f11+f14-f6-f16-f17-f12;
//........Store in "opposite" memory location..........
disteven[3*N+n] = f6;
disteven[6*N+n] = f12;
distodd[6*N+n] = f13;
disteven[8*N+n] = f16;
distodd[8*N+n] = f17;
//...................................................
}
}
__global__ void dvc_D3Q19_Velocity_BC_Z(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz, int outlet){
int n,N;
// distributions
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double uz;
N = Nx*Ny*Nz;
n = outlet + blockIdx.x*blockDim.x + threadIdx.x;
// Loop over the boundary - threadblocks delineated by start...finish
if ( n<N-Nx*Ny ){
// Read distributions from "opposite" memory convention
//........................................................................
f1 = distodd[n];
f3 = distodd[N+n];
f5 = distodd[2*N+n];
f7 = distodd[3*N+n];
f9 = distodd[4*N+n];
f11 = distodd[5*N+n];
f13 = distodd[6*N+n];
f15 = distodd[7*N+n];
f17 = distodd[8*N+n];
//........................................................................
f0 = disteven[n];
f2 = disteven[N+n];
f4 = disteven[2*N+n];
f6 = disteven[3*N+n];
f8 = disteven[4*N+n];
f10 = disteven[5*N+n];
f12 = disteven[6*N+n];
f14 = disteven[7*N+n];
f16 = disteven[8*N+n];
f18 = disteven[9*N+n];
//uz = -1.0 + (f0+f4+f3+f2+f1+f8+f7+f9+f10 + 2*(f6+f16+f17+f12+f13))/dout;
dout = (f0+f4+f3+f2+f1+f8+f7+f9+f10 + 2*(f6+f16+f17+f12+f13))/(1.0+uz);
f5 = f6 - 0.33333333333333338*dout* uz;
f15 = f16 - 0.16666666666666678*dout* uz;
f18 = f15 - f4 + f3-f16+f17-f8+f7-f9+f10;
f11 = (-dout*uz+f6+ f16+f17+f12+f13-f5-f15-f18+f2-f1-f13+f12+f8-f7-f9+f10)*0.5;
f14 = -dout*uz+f6+ f16+f17+f12+f13-f5-f15-f18-f11;
//........Store in "opposite" memory location..........
distodd[2*N+n] = f5;
distodd[5*N+n] = f11;
disteven[7*N+n] = f14;
distodd[7*N+n] = f15;
disteven[9*N+n] = f18;
//...................................................
}
}
extern "C" void PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N){
int GRID = count / 512 + 1;
@ -385,3 +495,17 @@ extern "C" void ComputePressureD3Q19(char *ID, double *disteven, double *distodd
int Nx, int Ny, int Nz){
dvc_ComputePressureD3Q19<<< NBLOCKS,NTHREADS >>>(ID, disteven, distodd, Pressure, Nx, Ny, Nz);
}
extern "C" void ScaLBL_D3Q19_Velocity_BC_z(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz){
int GRID = Nx*Ny / 512 + 1;
dvc_D3Q19_Velocity_BC_z<<<GRID,512>>>(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz);
}
extern "C" void ScaLBL_D3Q19_Velocity_BC_Z(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz, int outlet){
int GRID = Nx*Ny / 512 + 1;
dvc_D3Q19_Velocity_BC_Z<<<GRID,512>>>(double *disteven, double *distodd, double uz,
int Nx, int Ny, int Nz, int outlet);
}

10
tests/BasicSimulator.cpp
View File

@ -1103,7 +1103,7 @@ int main(int argc, char **argv)
// double vx_w_global,vy_w_global,vz_w_global; // global phase averaged velocity
// double vx_n_global,vy_n_global,vz_n_global; // global phase averaged velocity
double As_global;
double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0)
double dEs=0, dAwn=0, dAns=0; // Global surface energy (calculated by rank=0)
double pan_global,paw_global,pc_global; // local phase averaged pressure
DoubleArray van_global(3);
DoubleArray vaw_global(3);
@ -1439,12 +1439,12 @@ int main(int argc, char **argv)
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
//.......create and start timer............
double starttime,stoptime,cputime;
double starttime=0, stoptime=0, cputime=0;
sendtag = recvtag = 5;
FILE *TIMELOG;
FILE *FINALSTATE;
FILE *SPEED;
FILE *TIMELOG=NULL;
FILE *FINALSTATE=NULL;
FILE *SPEED=NULL;
if (rank==0){
TIMELOG= fopen("timelog.tcat","a+");
FINALSTATE= fopen("finalstate.tcat","a");

2
tests/ComponentLabel.cpp
View File

@ -422,7 +422,7 @@ int main(int argc, char **argv)
}
}
porosity /= (Nx*Ny*Nz*1.0);
printf("Media porosity is %f \n",porosity);
//printf("Media porosity is %f \n",porosity);
double beta=0.95;
int timestep=5;

234
tests/TestBlobIdentify.cpp
View File

@ -26,14 +26,17 @@ 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 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)));
inline double dist( const Point& p, double Lx, double Ly, double Lz ) const {
double x = fabs(p.x-center.x);
double y = fabs(p.y-center.y);
double z = fabs(p.z-center.z);
x = std::min(x,fabs(Lx-x));
y = std::min(y,fabs(Ly-y));
z = std::min(z,fabs(Lz-z));
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 +76,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 +138,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_ENABLE(1);
PROFILE_DISABLE_TRACE();
PROFILE_SYNCHRONIZE();
PROFILE_START("main");
#endif
if ( rank==0 ) {
printf("-----------------------------------------------------------\n");
printf("Testing Blob Identification \n");
@ -109,22 +164,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 +209,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 +218,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;
}

51
tests/lbpm_color_simulator.cpp
View File

@ -272,15 +272,18 @@ int main(int argc, char **argv)
printf("Parallel domain size = %i x %i x %i\n",nprocx,nprocy,nprocz);
if (BoundaryCondition==0) printf("Periodic boundary conditions will applied \n");
if (BoundaryCondition==1) printf("Pressure boundary conditions will be applied \n");
if (BoundaryCondition==2) printf("Velocity boundary conditions will be applied \n");
if (InitialCondition==0) printf("Initial conditions assigned from phase ID file \n");
if (InitialCondition==1) printf("Initial conditions asdsigned from restart file \n");
printf("********************************************************\n");
}
// Initialized domain and averaging framework for Two-Phase Flow
bool pBC;
bool pBC,velBC;
if (BoundaryCondition==1) pBC=true;
else pBC=false;
if (BoundaryCondition==2) velBC=true;
else velBC=false;
bool Restart;
if (InitialCondition==1) Restart=true;
else Restart=false;
@ -376,8 +379,8 @@ int main(int argc, char **argv)
int kstart,kfinish;
kstart = 1;
kfinish = Nz-1;
if (pBC && kproc==0) kstart = 4;
if (pBC && kproc==nprocz-1) kfinish = Nz-4;
if (BoundaryCondition > 0 && kproc==0) kstart = 4;
if (BoundaryCondition > 0 && kproc==nprocz-1) kfinish = Nz-4;
// Compute the pore volume
sum_local = 0.0;
@ -396,8 +399,8 @@ int main(int argc, char **argv)
porosity = pore_vol*iVol_global;
if (rank==0) printf("Media porosity = %f \n",porosity);
//.........................................................
// If pressure boundary conditions are applied remove solid
if (pBC && kproc == 0){
// If external boundary conditions are applied remove solid
if (BoundaryCondition > 0 && kproc == 0){
for (k=0; k<3; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
@ -408,7 +411,7 @@ int main(int argc, char **argv)
}
}
}
if (pBC && kproc == nprocz-1){
if (BoundaryCondition > 0 && kproc == nprocz-1){
for (k=Nz-3; k<Nz; k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
@ -541,20 +544,34 @@ int main(int argc, char **argv)
MPI_Barrier(MPI_COMM_WORLD);
//*************************************************************************
if (rank==0 && pBC){
if (rank==0 && BoundaryCondition==1){
printf("Setting inlet pressure = %f \n", din);
printf("Setting outlet pressure = %f \n", dout);
}
if (pBC && kproc == 0) {
if (BoundaryCondition==1 && kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (pBC && kproc == nprocz-1){
if (BoundaryCondition==1 && kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (rank==0 && BoundaryCondition==2){
printf("Setting inlet velocity = %f \n", din);
printf("Setting outlet velocity = %f \n", dout);
}
if (BoundaryCondition==2 && kproc == 0) {
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (BoundaryCondition==2 && kproc == nprocz-1){
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
ComputeVelocityD3Q19(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
@ -669,15 +686,25 @@ int main(int argc, char **argv)
DeviceBarrier();
if (pBC && kproc == 0) {
// Pressure boundary conditions
if (BoundaryCondition==1 && kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (pBC && kproc == nprocz-1){
if (BoundaryCondition==1 && kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
// Velocity boundary conditions
if (BoundaryCondition==2 && kproc == 0) {
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (BoundaryCondition==2 && kproc == nprocz-1){
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
//...................................................................................
MPI_Barrier(MPI_COMM_WORLD);

10
tests/testUtilities.cpp
View File

@ -10,6 +10,7 @@
#include <stdint.h>
#include "common/Utilities.h"
#include "common/StackTrace.h"
#include "common/UnitTest.h"
#include "common/MPI_Helpers.h"
@ -29,10 +30,13 @@
// Function to return the call stack
std::vector<std::string> get_call_stack()
{
std::vector<std::string> stack = Utilities::getCallStack();
std::vector<StackTrace::stack_info> stack = StackTrace::getCallStack();
std::vector<std::string> stack2(stack.size());
for (size_t i=0; i<stack.size(); i++)
stack2[i] = stack[i].print();
// Trick compiler to skip inline for this function with fake recursion
if ( stack.size() > 10000 ) { stack = get_call_stack(); }
return stack;
if ( stack.size() > 10000 ) { stack2 = get_call_stack(); }
return stack2;
}