From 63ef7d0fc699e5639e3aaeed9b7509dbff46da6e Mon Sep 17 00:00:00 2001 From: Mark Berrill Date: Tue, 25 Aug 2015 14:32:36 -0400 Subject: [PATCH] Finishing id mapping between timesteps, modifying TwoPhase to print the new global id. Still testing. --- README.titan | 4 +- analysis/analysis.cpp | 253 ++++-- analysis/analysis.h | 71 +- cmake/libraries.cmake | 11 +- common/Domain.h | 10 +- common/TwoPhase.cpp | 1392 ++++++++++++++++++++++++++++++++ common/TwoPhase.h | 1262 +---------------------------- common/pmmc.h | 16 +- pmmc/pmmc.cpp | 26 +- tests/CMakeLists.txt | 2 +- tests/TestBlobIdentify.cpp | 32 +- tests/lbpm_color_simulator.cpp | 8 +- tests/lbpm_color_simulator.h | 56 +- visit/avtLBMFileFormat.C | 18 +- visit/vtkHelpers.h | 2 +- 15 files changed, 1792 insertions(+), 1371 deletions(-) create mode 100644 common/TwoPhase.cpp diff --git a/README.titan b/README.titan index 0e365ed1..8d087700 100644 --- a/README.titan +++ b/README.titan @@ -24,8 +24,8 @@ cmake \ -D CMAKE_C_COMPILER:PATH=cc \ -D CMAKE_CXX_COMPILER:PATH=CC \ -D CMAKE_CXX_COMPILER:PATH=CC \ - -D CMAKE_C_FLAGS="-DCBUB" \ - -D CMAKE_CXX_FLAGS="-DCBUB" \ + -D CFLAGS="-DCBUB" \ + -D CXXFLAGS="-DCBUB" \ -D MPI_COMPILER:BOOL=TRUE \ -D MPIEXEC=aprun \ -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \ diff --git a/analysis/analysis.cpp b/analysis/analysis.cpp index 23e707bb..87a037aa 100644 --- a/analysis/analysis.cpp +++ b/analysis/analysis.cpp @@ -12,7 +12,7 @@ inline const TYPE* getPtr( const std::vector& x ) { return x.empty() ? NUL /****************************************************************** * Compute the blobs * ******************************************************************/ -int ComputeBlob( const Array& isPhase, IntArray& LocalBlobID, bool periodic, int start_id ) +int ComputeBlob( const Array& isPhase, BlobIDArray& LocalBlobID, bool periodic, int start_id ) { PROFILE_START("ComputeBlob",1); ASSERT(isPhase.size()==LocalBlobID.size()); @@ -53,7 +53,7 @@ int ComputeBlob( const Array& isPhase, IntArray& LocalBlobID, bool periodi std::vector neighbor_ids; neighbor_ids.reserve(N_neighbors); const bool *isPhasePtr = isPhase.get(); - int *LocalBlobIDPtr = LocalBlobID.get(); + BlobIDType *LocalBlobIDPtr = LocalBlobID.get(); for (int z=0; z& isPhase, IntArray& LocalBlobID, bool periodi * Compute the local blob ids * ******************************************************************/ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist, - double vF, double vS, IntArray& LocalBlobID, bool periodic ) + double vF, double vS, BlobIDArray& LocalBlobID, bool periodic ) { PROFILE_START("ComputeLocalBlobIDs"); ASSERT(SignDist.size()==Phase.size()); @@ -158,7 +158,7 @@ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist, PROFILE_STOP("ComputeLocalBlobIDs"); return nblobs; } -int ComputeLocalPhaseComponent(const IntArray &PhaseID, int VALUE, IntArray &ComponentLabel, bool periodic ) +int ComputeLocalPhaseComponent(const IntArray &PhaseID, int VALUE, BlobIDArray &ComponentLabel, bool periodic ) { PROFILE_START("ComputeLocalPhaseComponent"); size_t Nx = PhaseID.size(0); @@ -186,7 +186,7 @@ int ComputeLocalPhaseComponent(const IntArray &PhaseID, int VALUE, IntArray &Com /****************************************************************** * Reorder the global blob ids * ******************************************************************/ -static int ReorderBlobIDs2( IntArray& ID, int N_blobs, int ngx, int ngy, int ngz ) +static int ReorderBlobIDs2( BlobIDArray& ID, int N_blobs, int ngx, int ngy, int ngz ) { if ( N_blobs==0 ) return 0; @@ -233,7 +233,7 @@ static int ReorderBlobIDs2( IntArray& ID, int N_blobs, int ngx, int ngy, int ngz PROFILE_STOP("ReorderBlobIDs2",1); return N_blobs2; } -void ReorderBlobIDs( IntArray& ID ) +void ReorderBlobIDs( BlobIDArray& ID ) { PROFILE_START("ReorderBlobIDs"); int tmp = ID.max()+1; @@ -301,7 +301,7 @@ static bool updateLocalIds( const std::map& remote_map, return changed; } static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info, - int nblobs, IntArray& IDs ) + int nblobs, BlobIDArray& IDs ) { PROFILE_START("LocalToGlobalIDs",1); const int rank = rank_info.rank[1][1][1]; @@ -327,9 +327,9 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_ if ( IDs(i) >= 0 ) IDs(i) += offset; } - const Array LocalIDs = IDs; + const BlobIDArray LocalIDs = IDs; // Copy the ids and get the neighbors through the halos - fillHalo fillData(rank_info,nx,ny,nz,1,1,1,0,1,true,true,true); + fillHalo 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 neighbors; @@ -414,14 +414,14 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_ for (size_t k=ngz; k= 0 ) IDs(i,j,k) = final_map[id-offset]; } } } // Fill the ghosts - fillHalo fillData2(rank_info,nx,ny,nz,1,1,1,0,1,true,true,true); + fillHalo 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); @@ -434,7 +434,7 @@ static int LocalToGlobalIDs( int nx, int ny, int nz, const RankInfoStruct& rank_ } int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info, const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS, - IntArray& GlobalBlobID ) + BlobIDArray& GlobalBlobID ) { PROFILE_START("ComputeGlobalBlobIDs"); // First compute the local ids @@ -445,7 +445,7 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_inf return nglobal; } int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& rank_info, - const IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID ) + const IntArray &PhaseID, int VALUE, BlobIDArray &GlobalBlobID ) { PROFILE_START("ComputeGlobalPhaseComponent"); // First compute the local ids @@ -460,34 +460,48 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& r /****************************************************************** * Compute the mapping of blob ids between timesteps * ******************************************************************/ -void gatherSet( std::set& set ) +template inline MPI_Datatype getMPIType(); +template<> inline MPI_Datatype getMPIType() { return MPI_INT; } +template<> inline MPI_Datatype getMPIType() { + if ( sizeof(int64_t)==sizeof(long int) ) + return MPI_LONG; + else if ( sizeof(int64_t)==sizeof(double) ) + return MPI_DOUBLE; +} +template +void gatherSet( std::set& set ) { int nprocs; MPI_Comm_size(MPI_COMM_WORLD,&nprocs); - std::vector send_data(set.begin(),set.end()); + MPI_Datatype type = getMPIType(); + std::vector send_data(set.begin(),set.end()); int send_count = send_data.size(); std::vector 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 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, + std::vector recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]); + MPI_Allgatherv(getPtr(send_data),send_count,type, + getPtr(recv_data),getPtr(recv_count),getPtr(recv_disp),type, MPI_COMM_WORLD); for (size_t i=0; i >& src_map ) +template +void gatherSrcIDMap( std::map >& src_map ) { int nprocs; MPI_Comm_size(MPI_COMM_WORLD,&nprocs); - std::vector send_data; - for (std::map >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) { + MPI_Datatype type = getMPIType(); + std::vector send_data; + typename std::map >::const_iterator it; + for (it=src_map.begin(); it!=src_map.end(); ++it) { int id = it->first; - const std::set& src_ids = it->second; + const std::set& src_ids = it->second; send_data.push_back(id); send_data.push_back(src_ids.size()); - for (std::set::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2) + typename std::set::const_iterator it2; + for (it2=src_ids.begin(); it2!=src_ids.end(); ++it2) send_data.push_back(*it2); } int send_count = send_data.size(); @@ -495,45 +509,45 @@ void gatherSrcIDMap( std::map >& src_map ) MPI_Allgather(&send_count,1,MPI_INT,getPtr(recv_count),1,MPI_INT,MPI_COMM_WORLD); for (int i=1; i 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, + std::vector recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]); + MPI_Allgatherv(getPtr(send_data),send_count,type, + getPtr(recv_data),getPtr(recv_count),getPtr(recv_disp),type, 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& src_ids = src_map[id]; + std::set& src_ids = src_map[id]; for (int j=0; j >& src_map, - std::map >& dst_map, std::set& src, std::set& dst ) +void addSrcDstIDs( BlobIDType src_id, std::map >& src_map, + std::map >& dst_map, std::set& src, std::set& dst ) { src.insert(src_id); - const std::set& dst_ids = dst_map[src_id]; - for (std::set::const_iterator it=dst_ids.begin(); it!=dst_ids.end(); ++it) { + const std::set& dst_ids = dst_map[src_id]; + for (std::set::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 ) +ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2 ) { ASSERT(ID1.size()==ID2.size()); PROFILE_START("computeIDMap"); // Get a global list of all src/dst ids and the src map for each local blob - std::set src_set, dst_set; - std::map > src_map; // Map of the src ids for each dst id + std::set src_set, dst_set; + std::map > src_map; // Map of the src ids for each dst id for (size_t i=0; 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& src_ids = src_map[ID2(i)]; + std::set& src_ids = src_map[ID2(i)]; src_ids.insert(ID1(i)); } } @@ -542,12 +556,12 @@ ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 ) gatherSet( dst_set ); gatherSrcIDMap( src_map ); // Compute the dst id map - std::map > dst_map; // Map of the dst ids for each src id - for (std::map >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) { - int id = it->first; - const std::set& src_ids = it->second; - for (std::set::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2) { - std::set& dst_ids = dst_map[*it2]; + std::map > dst_map; // Map of the dst ids for each src id + for (std::map >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) { + BlobIDType id = it->first; + const std::set& src_ids = it->second; + for (std::set::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2) { + std::set& dst_ids = dst_map[*it2]; dst_ids.insert(id); } } @@ -555,53 +569,53 @@ ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 ) // Perform the mapping of ids ID_map_struct id_map; // Find new blobs - for (std::set::const_iterator it=dst_set.begin(); it!=dst_set.end(); ++it) { + for (std::set::const_iterator it=dst_set.begin(); it!=dst_set.end(); ++it) { if ( src_map.find(*it)==src_map.end() ) id_map.created.push_back(*it); } // Fine blobs that disappeared - for (std::set::const_iterator it=src_set.begin(); it!=src_set.end(); ++it) { + for (std::set::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 dst_list; + std::vector dst_list; dst_list.reserve(src_map.size()); - for (std::map >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) + for (std::map >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) dst_list.push_back(it->first); for (size_t i=0; i& src_ids = src_map[dst_id]; + const std::set& src_ids = src_map[dst_id]; if ( src_ids.size()==1 ) { int src_id = *src_ids.begin(); - const std::set& dst_ids = dst_map[src_id]; + const std::set& 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(src_id,dst_id)); + id_map.src_dst.push_back(std::pair(src_id,dst_id)); } } } // Handle merge/splits while ( !dst_map.empty() ) { // Get a lit of the src-dst ids - std::set src, dst; + std::set src, dst; addSrcDstIDs( dst_map.begin()->first, src_map, dst_map, src, dst ); - for (std::set::const_iterator it=src.begin(); it!=src.end(); ++it) + for (std::set::const_iterator it=src.begin(); it!=src.end(); ++it) dst_map.erase(*it); - for (std::set::const_iterator it=dst.begin(); it!=dst.end(); ++it) + for (std::set::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(dst.begin(),dst.end())) ); + id_map.split.push_back( BlobIDSplitStruct(*src.begin(),std::vector(dst.begin(),dst.end())) ); } else if ( dst.size()==1 ) { // Bubble merge - id_map.merge.push_back( BlobIDMergeStruct(std::vector(src.begin(),src.end()),*dst.begin()) ); + id_map.merge.push_back( BlobIDMergeStruct(std::vector(src.begin(),src.end()),*dst.begin()) ); } else { // Bubble split/merge id_map.merge_split.push_back( BlobIDMergeSplitStruct( - std::vector(src.begin(),src.end()), std::vector(dst.begin(),dst.end() ) ) ); + std::vector(src.begin(),src.end()), std::vector(dst.begin(),dst.end() ) ) ); } } ASSERT(src_map.empty()); @@ -611,4 +625,135 @@ ID_map_struct computeIDMap( const IntArray& ID1, const IntArray& ID2 ) } +/****************************************************************** +* Renumber the ids * +******************************************************************/ +void getNewIDs( ID_map_struct& map, BlobIDType& id_max, std::vector& new_ids ) +{ + new_ids.resize(0); + // Renumber the ids that map directly + for (size_t i=0; i(id1+1) ) + new_ids.resize(id1+1,-1); + new_ids[id1] = id2; + } + // Renumber the created blobs to create new ids + for (size_t i=0; i(id1+1) ) + new_ids.resize(id1+1,-1); + new_ids[id1] = id2; + } + // Renumber the blob splits to create new ids + for (size_t i=0; i(id1+1) ) + new_ids.resize(id1+1,-1); + new_ids[id1] = id2; + } + } + // Renumber the blob merges to create a new id + for (size_t i=0; i(id1+1) ) + new_ids.resize(id1+1,-1); + new_ids[id1] = id2; + } + // Renumber the blob merge/splits to create new ids + for (size_t i=0; i(id1+1) ) + new_ids.resize(id1+1,-1); + new_ids[id1] = id2; + } + } + +} +void renumberIDs( const std::vector& new_ids, BlobIDArray& IDs ) +{ + size_t N = IDs.length(); + BlobIDType* ids = IDs.get(); + for (size_t i=0; i=0 ) + ids[i] = new_ids[id]; + } +} + + +/****************************************************************** +* Write the id map for the given timestep * +******************************************************************/ +void writeIDMap( const ID_map_struct& map, long long int timestep, const std::string& filename ) +{ + int rank; + MPI_Comm_rank(MPI_COMM_WORLD,&rank); + if ( rank!=0 ) + return; + bool empty = map.created.empty() && map.destroyed.empty() && + map.split.empty() && map.merge.empty() && map.merge_split.empty(); + for (size_t i=0; i(map.created[i])); + for (size_t i=0; i(map.destroyed[i])); + for (size_t i=0; i(map.src_dst[i].first), + static_cast(map.src_dst[i].second)); + } + for (size_t i=0; i(map.split[i].first), + static_cast(map.split[i].second[0])); + for (size_t j=1; j(map.split[i].second[j])); + } + for (size_t i=0; i(map.merge[i].first[0])); + for (size_t j=1; j(map.merge[i].first[j])); + fprintf(fid,"-%lli",static_cast(map.merge[i].second)); + } + for (size_t i=0; i(map.merge_split[i].first[0])); + for (size_t j=1; j(map.merge_split[i].first[j])); + fprintf(fid,"-%lli",static_cast(map.merge_split[i].second[0])); + for (size_t j=1; j(map.merge_split[i].second[j])); + } + fprintf(fid,"\n"); + fclose(fid); +} diff --git a/analysis/analysis.h b/analysis/analysis.h index f0d562bc..cb90e196 100644 --- a/analysis/analysis.h +++ b/analysis/analysis.h @@ -9,6 +9,11 @@ #include +// Define types to use for blob ids +typedef int32_t BlobIDType; +typedef Array BlobIDArray; + + /*! * @brief Compute the blob * @details Compute the blob (F>vf|S>vs) starting from (i,j,k) - oil blob @@ -22,7 +27,7 @@ * @return Returns the number of blobs */ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist, - double vF, double vS, IntArray& LocalBlobID, bool periodic=true ); + double vF, double vS, BlobIDArray& LocalBlobID, bool periodic=true ); /*! * @brief Compute blob of an arbitrary phase @@ -33,8 +38,7 @@ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist, * @param[out] ComponentLabel * @param[in] periodic */ -int ComputeLocalPhaseComponent( const IntArray &PhaseID, int VALUE, IntArray &ComponentLabel, - bool periodic ); +int ComputeLocalPhaseComponent( const IntArray &PhaseID, int VALUE, IntArray &ComponentLabel, bool periodic ); /*! @@ -54,7 +58,7 @@ int ComputeLocalPhaseComponent( const IntArray &PhaseID, int VALUE, IntArray &Co */ int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info, const DoubleArray& Phase, const DoubleArray& SignDist, double vF, double vS, - IntArray& GlobalBlobID ); + BlobIDArray& GlobalBlobID ); /*! @@ -71,7 +75,7 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_inf * @return Return the number of components in the specified phase */ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& rank_info, - const IntArray &PhaseID, int VALUE, IntArray &GlobalBlobID ); + const IntArray &PhaseID, int VALUE, BlobIDArray &GlobalBlobID ); /*! @@ -83,19 +87,26 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& r * @param[in] nz Number of elements in the z-direction * @param[in/out] ID The ids of the blobs */ -void ReorderBlobIDs( IntArray& ID ); +void ReorderBlobIDs( BlobIDArray& ID ); -typedef std::pair > BlobIDSplitStruct; -typedef std::pair,int> BlobIDMergeStruct; -typedef std::pair,std::vector > BlobIDMergeSplitStruct; +typedef std::pair > BlobIDSplitStruct; +typedef std::pair,BlobIDType> BlobIDMergeStruct; +typedef std::pair,std::vector > BlobIDMergeSplitStruct; struct ID_map_struct { - std::vector created; // list of new blobs that were created - std::vector destroyed; // list of blobs that disappeared - std::vector > src_dst; // one-one mapping of blobs (first,second timestep id) + std::vector created; // list of new blobs that were created + std::vector destroyed; // list of blobs that disappeared + std::vector > src_dst; // one-one mapping of blobs (first,second timestep id) std::vector split; // list of blobs that split std::vector merge; // list of blobs that merged std::vector merge_split; // list of blobs that both merged and split + //! Empty constructor + ID_map_struct() {} + //! Create initial map from N blobs (ordered 1:N-1) + ID_map_struct( int N ) { + created.resize(N); + for (int i=0; i& new_ids ); + + +/*! + * @brief Update the blob ids based on mapping + * @details This functions computes the map of blob ids between iterations. + * Note: we also update the map to reflect the new ids + * @param[out] new_ids The newly renumbered blob ids (0:ids.max()) + * @param[in/out] IDs The blob ids to renumber + */ +void renumberIDs( const std::vector& new_id_list, BlobIDArray& IDs ); + + +/*! + * @brief Write the ID map + * @details This functions writes the id map fo an iteration. + * If no ids changed, then nothing will be written + * Note: only rank 0 writes, and the file is created on timestep 0. + * @param[in] map The timestep mapping for the ids + * @param[in] timestep The current timestep (timestep 0 creates the file) + * @param[in] filename The filename to write/append + */ +void writeIDMap( const ID_map_struct& map, long long int timestep, const std::string& filename ); + #endif diff --git a/cmake/libraries.cmake b/cmake/libraries.cmake index b843c2e6..9ed75e5a 100644 --- a/cmake/libraries.cmake +++ b/cmake/libraries.cmake @@ -221,10 +221,17 @@ MACRO( CONFIGURE_SYSTEM ) SET( SYSTEM_LDFLAGS ${SYSTEM_LDFLAGS} -rdynamic ) ENDIF() ENDIF() + # Try to add -fPIC + SET( CMAKE_REQUIRED_FLAGS "${CMAKE_CXX_FLAGS} ${COVERAGE_FLAGS} -fPIC" ) + CHECK_CXX_SOURCE_COMPILES( "int main() { return 0;}" fPIC ) + IF ( fPIC ) + SET( SYSTEM_LDFLAGS ${SYSTEM_LDFLAGS} -fPIC ) + SET( SYSTEM_LDFLAGS ${SYSTEM_LDFLAGS} -fPIC ) + ENDIF() IF ( USING_GCC ) SET( SYSTEM_LIBS ${SYSTEM_LIBS} -lgfortran ) - SET(CMAKE_C_FLAGS " ${CMAKE_C_FLAGS} -fPIC" ) - SET(CMAKE_CXX_FLAGS " ${CMAKE_CXX_FLAGS} -fPIC" ) + SET(CFLAGS_EXTRA " ${CFLAGS_EXTRA} -fPIC" ) + SET(CXXFLAGS_EXTRA " ${CXXFLAGS_EXTRA} -fPIC" ) ENDIF() ELSEIF( ${CMAKE_SYSTEM_NAME} STREQUAL "Darwin" ) # Max specific system libraries diff --git a/common/Domain.h b/common/Domain.h index 4683470a..bd60c709 100755 --- a/common/Domain.h +++ b/common/Domain.h @@ -16,7 +16,7 @@ #include "common/MPI_Helpers.h" #include "common/Communication.h" -int MAX_BLOB_COUNT=50; +static int MAX_BLOB_COUNT=50; using namespace std; @@ -211,7 +211,7 @@ static inline void fgetl( char * str, int num, FILE * stream ) if ( 0 ) {char *temp = (char *)&ptr; temp++;} } -Domain::~Domain(){ +inline Domain::~Domain(){ delete sendData_x; delete sendData_y; delete sendData_z; @@ -250,7 +250,7 @@ Domain::~Domain(){ delete recvData_YZ; } -void Domain::InitializeRanks() +inline void Domain::InitializeRanks() { // map the rank to the block index iproc = rank%nprocx; @@ -282,7 +282,7 @@ void Domain::InitializeRanks() } -void Domain::CommInit(MPI_Comm Communicator){ +inline void Domain::CommInit(MPI_Comm Communicator){ int i,j,k,n; int sendtag = 21; int recvtag = 21; @@ -696,7 +696,7 @@ inline void Domain::CommunicateMeshHalo(DoubleArray &Mesh) UnpackMeshData(recvList_YZ, recvCount_YZ ,recvData_YZ, MeshData); } -void Domain::BlobComm(MPI_Comm Communicator){ +inline void Domain::BlobComm(MPI_Comm Communicator){ //...................................................................................... int sendtag, recvtag; sendtag = recvtag = 51; diff --git a/common/TwoPhase.cpp b/common/TwoPhase.cpp new file mode 100644 index 00000000..476a5ccc --- /dev/null +++ b/common/TwoPhase.cpp @@ -0,0 +1,1392 @@ +#include "TwoPhase.h" + +#include "pmmc.h" +#include "Domain.h" +#include "Communication.h" +#include "analysis/analysis.h" + +#include "shared_ptr.h" +#include "common/Utilities.h" +#include "common/MPI_Helpers.h" +#include "IO/MeshDatabase.h" +#include "IO/Reader.h" +#include "IO/Writer.h" + +#define BLOB_AVG_COUNT 33 + +// Array access for averages defined by the following +#define VOL 0 +#define TRIMVOL 1 +#define PRS 2 +#define AWN 3 +#define AWS 4 +#define ANS 5 +#define LWNS 6 +#define JWN 7 +#define KWN 8 +#define CWNS 9 +#define KNWNS 10 +#define KGWNS 11 +#define VX 12 +#define VY 13 +#define VZ 14 +#define VSQ 15 +#define VWNX 16 +#define VWNY 17 +#define VWNZ 18 +#define VWNSX 19 +#define VWNSY 20 +#define VWNSZ 21 +#define GWNXX 22 +#define GWNYY 23 +#define GWNZZ 24 +#define GWNXY 25 +#define GWNXZ 26 +#define GWNYZ 27 +#define TRAWN 28 +#define TRJWN 29 +#define CMX 30 +#define CMY 31 +#define CMZ 32 + + +// Constructor +TwoPhase::TwoPhase(Domain &dm) : Dm(dm) +{ + Nx=dm.Nx; Ny=dm.Ny; Nz=dm.Nz; + Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz*1.0; + + // Global arrays + PhaseID.resize(Nx,Ny,Nz); + Label_WP.resize(Nx,Ny,Nz); + Label_NWP.resize(Nx,Ny,Nz); + Label_NWP.resize(Nx,Ny,Nz); + SDn.resize(Nx,Ny,Nz); + SDs.resize(Nx,Ny,Nz); + Phase.resize(Nx,Ny,Nz); + Press.resize(Nx,Ny,Nz); + dPdt.resize(Nx,Ny,Nz); + MeanCurvature.resize(Nx,Ny,Nz); + GaussCurvature.resize(Nx,Ny,Nz); + SDs_x.resize(Nx,Ny,Nz); // Gradient of the signed distance + SDs_y.resize(Nx,Ny,Nz); + SDs_z.resize(Nx,Ny,Nz); + SDn_x.resize(Nx,Ny,Nz); // Gradient of the signed distance + SDn_y.resize(Nx,Ny,Nz); + SDn_z.resize(Nx,Ny,Nz); + DelPhi.resize(Nx,Ny,Nz); + Phase_tplus.resize(Nx,Ny,Nz); + Phase_tminus.resize(Nx,Ny,Nz); + Vel_x.resize(Nx,Ny,Nz); // Gradient of the phase indicator field + Vel_y.resize(Nx,Ny,Nz); + Vel_z.resize(Nx,Ny,Nz); + //......................................... + // Allocate cube storage space + CubeValues.resize(2,2,2); + nw_tris.resize(3,20); + ns_tris.resize(3,20); + ws_tris.resize(3,20); + nws_seg.resize(2,20); + local_sol_tris.resize(3,18); + nw_pts=DTMutableList(20); + ns_pts=DTMutableList(20); + ws_pts=DTMutableList(20); + nws_pts=DTMutableList(20); + local_nws_pts=DTMutableList(20); + local_sol_pts=DTMutableList(20); + tmp=DTMutableList(20); + //......................................... + Values.resize(20); + DistanceValues.resize(20); + KGwns_values.resize(20); + KNwns_values.resize(20); + InterfaceSpeed.resize(20); + NormalVector.resize(60); + //......................................... + van.resize(3); + vaw.resize(3); + vawn.resize(3); + vawns.resize(3); + Gwn.resize(6); + Gns.resize(6); + Gws.resize(6); + van_global.resize(3); + vaw_global.resize(3); + vawn_global.resize(3); + vawns_global.resize(3); + Gwn_global.resize(6); + Gns_global.resize(6); + Gws_global.resize(6); + //......................................... + if (Dm.rank==0) +{ + TIMELOG = fopen("timelog.tcat","a+"); + if (fseek(TIMELOG,0,SEEK_SET) == fseek(TIMELOG,0,SEEK_CUR)) +{ + // If timelog is empty, write a short header to list the averages + //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); + fprintf(TIMELOG,"time dEs "); // Timestep, Change in Surface Energy + fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns cwns KNwns KGwns "); // Scalar averages + fprintf(TIMELOG,"vawx vawy vawz vanx vany vanz "); // Velocity averages + fprintf(TIMELOG,"vawnx vawny vawnz vawnsx vawnsy vawnsz "); + fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors + fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz "); + fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz "); + fprintf(TIMELOG,"trawn trJwn trRwn\n"); // trimmed curvature for wn surface + //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); + } + + NWPLOG = fopen("components.NWP.tcat","a+"); + fprintf(NWPLOG,"time label vol pn awn ans Jwn Kwn lwns cwns "); + fprintf(NWPLOG,"vx vy vz vwnx vwny vwnz vwnsx vwnsy vwnsz vsq "); + fprintf(NWPLOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz Cx Cy Cz trawn trJwn\n"); + + WPLOG = fopen("components.WP.tcat","a+"); + fprintf(WPLOG,"time label vol pw awn ans Jwn Kwn lwns cwns "); + fprintf(WPLOG,"vx vy vz vwnx vwny vwnz vwnsx vwnsy vwnsz vsq "); + fprintf(WPLOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz Cx Cy Cz trawn trJwn\n"); + } +} + + +// Destructor +TwoPhase::~TwoPhase() +{ +} + + +void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, DoubleArray &DistData) +{ +/* double factor,temp,value; + factor=0.5/Beta; + for (int n=0; n 0.999 ) DistData[n] = 4.0; + else if (value < -0.999 ) DistData[n] = -4.0; + else DistData[n] = factor*log((1.0+value)/(1.0-value)); + if (DistData[n] > 1.0) DistData[n] = 1.0; + if (DistData[n] < -1.0) DistData[n] = -1.0; + } + // Initialize to -1,1 (segmentation) + for (int k=0; k 1.0) DistData(i,j,k) = 1.0; + else if (temp < -1.0) DistData(i,j,k) = -1.0; + else DistData(i,j,k) = temp; + } + } + } + + SSO(DistData,Dm.id,Dm,10); +*/ + for (int k=0; k 0 && Dm.kproc == 0) kmin=4; + if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4; + + for (k=kmin; k 0 ) +{ + nwp_volume += 0.125; + // velocity + van(0) += 0.125*Vel_x(n); + van(1) += 0.125*Vel_y(n); + van(2) += 0.125*Vel_z(n); + // volume the excludes the interfacial region + if (DelPhi(n) < 1e-4) +{ + vol_n += 0.125; + // pressure + pan += 0.125*Press(n); + + } + } + else{ + wp_volume += 0.125; + // velocity + vaw(0) += 0.125*Vel_x(n); + vaw(1) += 0.125*Vel_y(n); + vaw(2) += 0.125*Vel_z(n); + if (DelPhi(n) < 1e-4) +{ + // volume the excludes the interfacial region + vol_w += 0.125; + // pressure + paw += 0.125*Press(n); + + } + } + } + } + + //........................................................................... + // Construct the interfaces and common curve + pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue, + nw_pts, nw_tris, Values, ns_pts, ns_tris, ws_pts, ws_tris, + local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris, + n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris, + n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg, + i, j, k, Nx, Ny, Nz); + + // wn interface averages + if (n_nw_pts > 0) +{ + awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris); + Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); + Kwn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); + + // Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels) + pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues, + i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn); + + pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues, + i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn); + + // Compute the normal speed of the interface + pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris, + NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris); + } + // wns common curve averages + if (n_local_nws_pts > 0) +{ + efawns += pmmc_CubeContactAngle(CubeValues,Values,SDn_x,SDn_y,SDn_z,SDs_x,SDs_y,SDs_z, + local_nws_pts,i,j,k,n_local_nws_pts); + + pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z, + local_nws_pts,i,j,k,n_local_nws_pts); + + pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y, + SDs_z, KNwns_values, KGwns_values, KNwns, KGwns, + nws_pts, n_nws_pts, i, j, k); + + lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts); + } + + // Solid interface averagees + if (n_local_sol_tris > 0) +{ + As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris); + + // Compute the surface orientation and the interfacial area + ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris); + aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris); + } + //........................................................................... + } + } + } +} +void TwoPhase::AssignComponentLabels() +{ + int LabelNWP=1; + //int LabelWP=2; + // NOTE: labeling the wetting phase components is tricky! One sandstone media had over 800,000 components + //NumberComponents_WP = ComputeGlobalPhaseComponent(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,PhaseID,LabelWP,Label_WP); + // treat all wetting phase is connected + NumberComponents_WP=1; + for (int k=0; k 0 && Dm.kproc == 0) kmin=4; + if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4; + + for (k=kmin; k 0.0 && !(LabelNWP < 0) ) +{ + // volume + ComponentAverages_NWP(VOL,LabelNWP) += 0.125; + // velocity + ComponentAverages_NWP(VX,LabelNWP) += 0.125*Vel_x(n); + ComponentAverages_NWP(VY,LabelNWP) += 0.125*Vel_y(n); + ComponentAverages_NWP(VZ,LabelNWP) += 0.125*Vel_z(n); + // center of mass + ComponentAverages_NWP(CMX,LabelNWP) += 0.125*(i+cube[p][0]+Dm.iproc*Nx); + ComponentAverages_NWP(CMY,LabelNWP) += 0.125*(j+cube[p][1]+Dm.jproc*Ny); + ComponentAverages_NWP(CMZ,LabelNWP) += 0.125*(k+cube[p][2]+Dm.kproc*Nz); + + // twice the kinetic energy + ComponentAverages_NWP(VSQ,LabelNWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n)); + + // volume the for pressure averaging excludes the interfacial region + if (DelPhi(n) < 1e-4 ) +{ + ComponentAverages_NWP(TRIMVOL,LabelNWP) += 0.125; + ComponentAverages_NWP(PRS,LabelNWP) += 0.125*Press(n); + } + } + else if (!(LabelWP < 0)) +{ + ComponentAverages_WP(VOL,LabelWP) += 0.125; + // velocity + ComponentAverages_WP(VX,LabelWP) += 0.125*Vel_x(n); + ComponentAverages_WP(VY,LabelWP)+= 0.125*Vel_y(n); + ComponentAverages_WP(VZ,LabelWP) += 0.125*Vel_z(n); + // Center of mass + ComponentAverages_WP(CMX,LabelWP) += 0.125*(i+cube[p][0]+Dm.iproc*Nx); + ComponentAverages_WP(CMY,LabelWP) += 0.125*(j+cube[p][1]+Dm.jproc*Ny); + ComponentAverages_WP(CMZ,LabelWP) += 0.125*(k+cube[p][2]+Dm.kproc*Nz); + // twice the kinetic energy + ComponentAverages_WP(VSQ,LabelWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n)); + + // volume the for pressure averaging excludes the interfacial region + if (DelPhi(n) < 1e-4) +{ + ComponentAverages_WP(TRIMVOL,LabelWP) += 0.125; + ComponentAverages_WP(PRS,LabelWP) += 0.125*Press(n); + } + } + } + } + //........................................................................... + // Construct the interfaces and common curve + pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue, + nw_pts, nw_tris, Values, ns_pts, ns_tris, ws_pts, ws_tris, + local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris, + n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris, + n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg, + i, j, k, Nx, Ny, Nz); + + //........................................................................... + // wn interface averages + if (n_nw_pts>0 && LabelNWP >=0 && LabelWP >=0 ) +{ + // Mean curvature + TempLocal = pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); + ComponentAverages_WP(JWN,LabelWP) += TempLocal; + ComponentAverages_NWP(JWN,LabelNWP) += TempLocal; + + // Trimmed Mean curvature + pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues, + i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn); + ComponentAverages_WP(TRAWN,LabelWP) += trawn; + ComponentAverages_WP(TRJWN,LabelWP) += trJwn; + ComponentAverages_NWP(TRAWN,LabelNWP) += trawn; + ComponentAverages_NWP(TRJWN,LabelNWP) += trJwn; + + // Gaussian curvature + TempLocal = pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); + ComponentAverages_WP(KWN,LabelWP) += TempLocal; + ComponentAverages_NWP(KWN,LabelNWP) += TempLocal; + + // Compute the normal speed of the interface + pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris, + NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris); + ComponentAverages_WP(VWNX,LabelWP) += vawn(0); + ComponentAverages_WP(VWNY,LabelWP) += vawn(1); + ComponentAverages_WP(VWNZ,LabelWP) += vawn(2); + ComponentAverages_NWP(VWNX,LabelNWP) += vawn(0); + ComponentAverages_NWP(VWNY,LabelNWP) += vawn(1); + ComponentAverages_NWP(VWNZ,LabelNWP) += vawn(2); + + // Interfacial Area + TempLocal = pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris); + ComponentAverages_WP(AWN,LabelWP) += TempLocal; + ComponentAverages_NWP(AWN,LabelNWP) += TempLocal; + + ComponentAverages_WP(GWNXX,LabelWP) += Gwn(0); + ComponentAverages_WP(GWNYY,LabelWP) += Gwn(1); + ComponentAverages_WP(GWNZZ,LabelWP) += Gwn(2); + ComponentAverages_WP(GWNXY,LabelWP) += Gwn(3); + ComponentAverages_WP(GWNXZ,LabelWP) += Gwn(4); + ComponentAverages_WP(GWNYZ,LabelWP) += Gwn(5); + + ComponentAverages_NWP(GWNXX,LabelNWP) += Gwn(0); + ComponentAverages_NWP(GWNYY,LabelNWP) += Gwn(1); + ComponentAverages_NWP(GWNZZ,LabelNWP) += Gwn(2); + ComponentAverages_NWP(GWNXY,LabelNWP) += Gwn(3); + ComponentAverages_NWP(GWNXZ,LabelNWP) += Gwn(4); + ComponentAverages_NWP(GWNYZ,LabelNWP) += Gwn(5); + + } + //........................................................................... + // Common curve averages + if (n_local_nws_pts > 0 && LabelNWP >=0 && LabelWP >=0) +{ + // Contact angle + TempLocal = pmmc_CubeContactAngle(CubeValues,Values,SDn_x,SDn_y,SDn_z,SDs_x,SDs_y,SDs_z, + local_nws_pts,i,j,k,n_local_nws_pts); + ComponentAverages_WP(CWNS,LabelWP) += TempLocal; + ComponentAverages_NWP(CWNS,LabelNWP) += TempLocal; + + // Kinematic velocity of the common curve + pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z, + local_nws_pts,i,j,k,n_local_nws_pts); + ComponentAverages_WP(VWNSX,LabelWP) += vawns(0); + ComponentAverages_WP(VWNSY,LabelWP) += vawns(1); + ComponentAverages_WP(VWNSZ,LabelWP) += vawns(2); + ComponentAverages_NWP(VWNSX,LabelNWP) += vawns(0); + ComponentAverages_NWP(VWNSY,LabelNWP) += vawns(1); + ComponentAverages_NWP(VWNSZ,LabelNWP) += vawns(2); + + // Curvature of the common curve + pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y, + SDs_z, KNwns_values, KGwns_values, KNwns, KGwns, + nws_pts, n_nws_pts, i, j, k); + ComponentAverages_WP(KNWNS,LabelWP) += KNwns; + ComponentAverages_WP(KGWNS,LabelWP) += KGwns; + ComponentAverages_NWP(KNWNS,LabelNWP) += KNwns; + ComponentAverages_NWP(KGWNS,LabelNWP) += KGwns; + + // Length of the common curve + TempLocal = pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts); + ComponentAverages_NWP(LWNS,LabelNWP) += TempLocal; + } + //........................................................................... + // Solid interface averages + if (n_local_sol_pts > 0 && LabelWP >=0) +{ + As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris); + // Compute the surface orientation and the interfacial area + + TempLocal = pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris); + ComponentAverages_WP(AWS,LabelWP) += TempLocal; + } + if (n_ns_pts > 0 && LabelNWP >=0 ) +{ + TempLocal = pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris); + ComponentAverages_NWP(ANS,LabelNWP) += TempLocal; + } + //........................................................................... + + } + } + } + + if (Dm.rank==0) +{ + printf("Component averages computed locally -- reducing result... \n"); + } + // Globally reduce the non-wetting phase averages + RecvBuffer.resize(BLOB_AVG_COUNT,NumberComponents_NWP); + +/* for (int b=0; b 0.0) +{ + double Vn,pn,awn,ans,Jwn,Kwn,lwns,cwns,vsq; + + Vn = ComponentAverages_NWP(VOL,b); + awn = ComponentAverages_NWP(AWN,b); + ans = ComponentAverages_NWP(ANS,b); + van(0) = ComponentAverages_NWP(VX,b)/Vn; + van(1) = ComponentAverages_NWP(VY,b)/Vn; + van(2) = ComponentAverages_NWP(VZ,b)/Vn; + vsq = ComponentAverages_NWP(VSQ,b)/Vn; + NULL_USE(ans); + + if (ComponentAverages_NWP(TRIMVOL,b) > 0.0) +{ + pn = ComponentAverages_NWP(PRS,b)/ComponentAverages_NWP(TRIMVOL,b); + } + else pn = 0.0; + + if (awn != 0.0) +{ + Jwn = ComponentAverages_NWP(JWN,b)/awn; + Kwn = ComponentAverages_NWP(KWN,b)/awn; + vawn(0) = ComponentAverages_NWP(VWNSX,b)/awn; + vawn(1) = ComponentAverages_NWP(VWNSY,b)/awn; + vawn(2) = ComponentAverages_NWP(VWNSZ,b)/awn; + Gwn(0) = ComponentAverages_NWP(GWNXX,b)/awn; + Gwn(1) = ComponentAverages_NWP(GWNYY,b)/awn; + Gwn(2) = ComponentAverages_NWP(GWNZZ,b)/awn; + Gwn(3) = ComponentAverages_NWP(GWNXY,b)/awn; + Gwn(4) = ComponentAverages_NWP(GWNXZ,b)/awn; + Gwn(5) = ComponentAverages_NWP(GWNYZ,b)/awn; + } + else Jwn=Kwn=0.0; + + trawn = ComponentAverages_NWP(TRAWN,b); + trJwn = ComponentAverages_NWP(TRJWN,b); + if (trawn > 0.0) trJwn /= trawn; + + lwns = ComponentAverages_NWP(LWNS,b); + if (lwns != 0.0) +{ + cwns = ComponentAverages_NWP(CWNS,b)/lwns; + vawns(0) = ComponentAverages_NWP(VWNSX,b)/lwns; + vawns(1) = ComponentAverages_NWP(VWNSY,b)/lwns; + vawns(2) = ComponentAverages_NWP(VWNSZ,b)/lwns; + } + else cwns=0.0; + + ComponentAverages_NWP(PRS,b) = pn; + ComponentAverages_NWP(VX,b) = van(0); + ComponentAverages_NWP(VY,b) = van(1); + ComponentAverages_NWP(VZ,b) = van(2); + ComponentAverages_NWP(VSQ,b) = vsq; + + ComponentAverages_NWP(JWN,b) = Jwn; + ComponentAverages_NWP(KWN,b) = Kwn; + ComponentAverages_NWP(VWNX,b) = vawn(0); + ComponentAverages_NWP(VWNY,b) = vawn(1); + ComponentAverages_NWP(VWNZ,b) = vawn(2); + + ComponentAverages_NWP(GWNXX,b) = Gwn(0); + ComponentAverages_NWP(GWNYY,b) = Gwn(1); + ComponentAverages_NWP(GWNZZ,b) = Gwn(2); + ComponentAverages_NWP(GWNXY,b) = Gwn(3); + ComponentAverages_NWP(GWNXZ,b) = Gwn(4); + ComponentAverages_NWP(GWNYZ,b) = Gwn(5); + + ComponentAverages_NWP(CWNS,b) = cwns; + ComponentAverages_NWP(VWNSX,b) = vawns(0); + ComponentAverages_NWP(VWNSY,b) = vawns(1); + ComponentAverages_NWP(VWNSZ,b) = vawns(2); + + ComponentAverages_NWP(CMX,b) /= Vn; + ComponentAverages_NWP(CMY,b) /= Vn; + ComponentAverages_NWP(CMZ,b) /= Vn; + + ComponentAverages_NWP(TRJWN,b) = trJwn; + + } + } + + // reduce the wetting phase averages + for (int b=0; b 0.0) +{ + double Vw,pw,awn,ans,Jwn,Kwn,lwns,cwns,vsq; + Vw = ComponentAverages_WP(VOL,b); + awn = ComponentAverages_WP(AWN,b); + ans = ComponentAverages_WP(ANS,b); + vaw(0) = ComponentAverages_WP(VX,b)/Vw; + vaw(1) = ComponentAverages_WP(VY,b)/Vw; + vaw(2) = ComponentAverages_WP(VZ,b)/Vw; + vsq = ComponentAverages_WP(VSQ,b)/Vw; + NULL_USE(ans); + + if (ComponentAverages_WP(TRIMVOL,b) > 0.0) +{ + pw = ComponentAverages_WP(PRS,b)/ComponentAverages_WP(TRIMVOL,b); + } + else pw = 0.0; + + if (awn != 0.0) +{ + Jwn = ComponentAverages_WP(JWN,b)/awn; + Kwn = ComponentAverages_WP(KWN,b)/awn; + vawn(0) = ComponentAverages_WP(VWNSX,b)/awn; + vawn(1) = ComponentAverages_WP(VWNSY,b)/awn; + vawn(2) = ComponentAverages_WP(VWNSZ,b)/awn; + Gwn(0) = ComponentAverages_WP(GWNXX,b)/awn; + Gwn(1) = ComponentAverages_WP(GWNYY,b)/awn; + Gwn(2) = ComponentAverages_WP(GWNZZ,b)/awn; + Gwn(3) = ComponentAverages_WP(GWNXY,b)/awn; + Gwn(4) = ComponentAverages_WP(GWNXZ,b)/awn; + Gwn(5) = ComponentAverages_WP(GWNYZ,b)/awn; + } + else Jwn=Kwn=0.0; + + trawn = ComponentAverages_WP(TRAWN,b); + trJwn = ComponentAverages_WP(TRJWN,b); + if (trawn > 0.0) trJwn /= trawn; + + lwns = ComponentAverages_WP(LWNS,b); + if (lwns != 0.0) +{ + cwns = ComponentAverages_WP(CWNS,b)/lwns; + vawns(0) = ComponentAverages_WP(VWNSX,b)/lwns; + vawns(1) = ComponentAverages_WP(VWNSY,b)/lwns; + vawns(2) = ComponentAverages_WP(VWNSZ,b)/lwns; + } + else cwns=0.0; + + ComponentAverages_WP(PRS,b) = pw; + ComponentAverages_WP(VX,b) = vaw(0); + ComponentAverages_WP(VY,b) = vaw(1); + ComponentAverages_WP(VZ,b) = vaw(2); + ComponentAverages_WP(VSQ,b) = vsq; + + ComponentAverages_WP(JWN,b) = Jwn; + ComponentAverages_WP(KWN,b) = Kwn; + ComponentAverages_WP(VWNX,b) = vawn(0); + ComponentAverages_WP(VWNY,b) = vawn(1); + ComponentAverages_WP(VWNZ,b) = vawn(2); + + ComponentAverages_WP(GWNXX,b) = Gwn(0); + ComponentAverages_WP(GWNYY,b) = Gwn(1); + ComponentAverages_WP(GWNZZ,b) = Gwn(2); + ComponentAverages_WP(GWNXY,b) = Gwn(3); + ComponentAverages_WP(GWNXZ,b) = Gwn(4); + ComponentAverages_WP(GWNYZ,b) = Gwn(5); + + ComponentAverages_WP(CWNS,b) = cwns; + ComponentAverages_WP(VWNSX,b) = vawns(0); + ComponentAverages_WP(VWNSY,b) = vawns(1); + ComponentAverages_WP(VWNSZ,b) = vawns(2); + + ComponentAverages_WP(TRJWN,b) = trJwn; + } + } +} + +void TwoPhase::WriteSurfaces(int logcount) +{ + int i,j,k; + int ncubes=(Nx-1)*(Ny-1)*(Nz-1); + Point P,A,B,C; + + std::shared_ptr wn_mesh( new IO::TriList() ); + wn_mesh->A.reserve(8*ncubes); + wn_mesh->B.reserve(8*ncubes); + wn_mesh->C.reserve(8*ncubes); + + std::shared_ptr ns_mesh( new IO::TriList() ); + ns_mesh->A.reserve(8*ncubes); + ns_mesh->B.reserve(8*ncubes); + ns_mesh->C.reserve(8*ncubes); + + std::shared_ptr ws_mesh( new IO::TriList() ); + ws_mesh->A.reserve(8*ncubes); + ws_mesh->B.reserve(8*ncubes); + ws_mesh->C.reserve(8*ncubes); + + for (k=1; kA.push_back(A); + wn_mesh->B.push_back(B); + wn_mesh->C.push_back(C); + } + for (int r=0;rA.push_back(A); + ws_mesh->B.push_back(B); + ws_mesh->C.push_back(C); + } + for (int r=0;rA.push_back(A); + ns_mesh->B.push_back(B); + ns_mesh->C.push_back(C); + } + } + } + } + + std::vector meshData(4); + meshData[0].meshName = "wn-tris"; + meshData[0].mesh = wn_mesh; + meshData[1].meshName = "ws-tris"; + meshData[1].mesh = ws_mesh; + meshData[2].meshName = "ns-tris"; + meshData[2].mesh = ns_mesh; + IO::writeData( logcount, meshData, 2); + +} + +void TwoPhase::Reduce() +{ + int i; + double iVol_global=1.0/Volume; + //........................................................................... + MPI_Barrier(Dm.Comm); + MPI_Allreduce(&nwp_volume,&nwp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&wp_volume,&wp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&awn,&awn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&ans,&ans_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&aws,&aws_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&lwns,&lwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&As,&As_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&Jwn,&Jwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&Kwn,&Kwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&KGwns,&KGwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&KNwns,&KNwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + // Phase averages + MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&paw,&paw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&pan,&pan_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&vaw(0),&vaw_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&van(0),&van_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&vawn(0),&vawn_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&vawns(0),&vawns_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&Gwn(0),&Gwn_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&Gns(0),&Gns_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&Gws(0),&Gws_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&trawn,&trawn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&trJwn,&trJwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Allreduce(&trRwn,&trRwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); + MPI_Barrier(Dm.Comm); + + // Normalize the phase averages + // (density of both components = 1.0) + if (vol_w_global > 0.0) +{ + paw_global = paw_global / vol_w_global; + } + if (wp_volume_global > 0.0) +{ + vaw_global(0) = vaw_global(0) / wp_volume_global; + vaw_global(1) = vaw_global(1) / wp_volume_global; + vaw_global(2) = vaw_global(2) / wp_volume_global; + + } + if (vol_n_global > 0.0) +{ + pan_global = pan_global / vol_n_global; + } + if (nwp_volume_global > 0.0) +{ + van_global(0) = van_global(0) / nwp_volume_global; + van_global(1) = van_global(1) / nwp_volume_global; + van_global(2) = van_global(2) / nwp_volume_global; + } + // Normalize surface averages by the interfacial area + if (awn_global > 0.0) +{ + Jwn_global /= awn_global; + Kwn_global /= awn_global; + for (i=0; i<3; i++) vawn_global(i) /= awn_global; + for (i=0; i<6; i++) Gwn_global(i) /= awn_global; + } + if (lwns_global > 0.0) +{ + efawns_global /= lwns_global; + KNwns_global /= lwns_global; + KGwns_global /= lwns_global; + for (i=0; i<3; i++) vawns_global(i) /= lwns_global; + } + if (trawn_global > 0.0) +{ + trJwn_global /= trawn_global; + trRwn_global /= trawn_global; + trRwn_global = 2.0*fabs(trRwn_global); + trJwn_global = fabs(trJwn_global); + } + + if (ans_global > 0.0) for (i=0; i<6; i++) Gns_global(i) /= ans_global; + if (aws_global > 0.0) for (i=0; i<6; i++) Gws_global(i) /= aws_global; + + //sat_w = 1.0 - nwp_volume_global*iVol_global/porosity; + sat_w = 1.0 - nwp_volume_global/(nwp_volume_global+wp_volume_global); + // Compute the specific interfacial areas and common line length (dimensionless per unit volume) + awn_global = awn_global*iVol_global; + ans_global = ans_global*iVol_global; + aws_global = aws_global*iVol_global; + dEs = dEs*iVol_global; + lwns_global = lwns_global*iVol_global; +} + +void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) +{ + awn_global *= D; + ans_global *= D; + ans_global *= D; + lwns_global *= D*D; +} + +void TwoPhase::PrintAll(int timestep) +{ + if (Dm.rank==0) +{ + fprintf(TIMELOG,"%i %.5g ",timestep-5,dEs); // change in surface energy + fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure + fprintf(TIMELOG,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas + fprintf(TIMELOG,"%.5g %.5g ",Jwn_global, Kwn_global); // curvature of wn interface + fprintf(TIMELOG,"%.5g ",lwns_global); // common curve length + fprintf(TIMELOG,"%.5g ",efawns_global); // average contact angle + fprintf(TIMELOG,"%.5g %.5g ",KNwns_global, KGwns_global); // curvature of wn interface + fprintf(TIMELOG,"%.5g %.5g %.5g ",vaw_global(0),vaw_global(1),vaw_global(2)); // average velocity of w phase + fprintf(TIMELOG,"%.5g %.5g %.5g ",van_global(0),van_global(1),van_global(2)); // average velocity of n phase + fprintf(TIMELOG,"%.5g %.5g %.5g ",vawn_global(0),vawn_global(1),vawn_global(2)); // velocity of wn interface + fprintf(TIMELOG,"%.5g %.5g %.5g ",vawns_global(0),vawns_global(1),vawns_global(2)); // velocity of wn interface + fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", + Gwn_global(0),Gwn_global(1),Gwn_global(2),Gwn_global(3),Gwn_global(4),Gwn_global(5)); // orientation of wn interface + fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", + Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5)); // orientation of ns interface + fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", + Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface + fprintf(TIMELOG,"%.5g %.5g %.5g\n",trawn_global, trJwn_global, trRwn_global); // Trimmed curvature + fflush(TIMELOG); + } +} + +void TwoPhase::PrintComponents(int timestep) +{ + if (Dm.rank==0) +{ + printf("PRINT %i COMPONENT AVEREAGES: time = %i \n",(int)ComponentAverages_NWP.size(1),timestep); + for (int b=0; b 0.0) +{ + fprintf(NWPLOG,"%i ",timestep-5); + if ( Label_NWP_map.empty() ) { + // print index + fprintf(NWPLOG,"%i ",b); + } else { + // print final global id + fprintf(NWPLOG,"%i ",Label_NWP_map[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 ",ComponentAverages_NWP(GWNYZ,b)); + fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CMX,b)); + fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CMY,b)); + fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CMZ,b)); + fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(TRAWN,b)); + fprintf(NWPLOG,"%.5g\n",ComponentAverages_NWP(TRJWN,b)); + } + } + fflush(NWPLOG); + + for (int b=0; 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 ",ComponentAverages_WP(GWNYZ,b)); + fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CMX,b)); + fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CMY,b)); + fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CMZ,b)); + fprintf(WPLOG,"%.5g ",ComponentAverages_WP(TRAWN,b)); + fprintf(WPLOG,"%.5g\n",ComponentAverages_WP(TRJWN,b)); + + } + } + fflush(WPLOG); + } +} + +inline int TwoPhase::GetCubeLabel(int i, int j, int k, IntArray &BlobLabel) +{ + int label; + label=BlobLabel(i,j,k); + label=max(label,BlobLabel(i+1,j,k)); + label=max(label,BlobLabel(i,j+1,k)); + label=max(label,BlobLabel(i+1,j+1,k)); + label=max(label,BlobLabel(i,j,k+1)); + label=max(label,BlobLabel(i+1,j,k+1)); + label=max(label,BlobLabel(i,j+1,k+1)); + label=max(label,BlobLabel(i+1,j+1,k+1)); + + return label; +} + +void TwoPhase::SortBlobs() +{ + //printf("Sorting the blobs based on volume \n"); + //printf("-----------------------------------------------\n"); + int TempLabel,a,aa,bb,i,j,k,idx; + double TempValue; + //....................................................................... + // Sort NWP components by volume + //....................................................................... + IntArray OldLabel(NumberComponents_NWP); + for (a=0; a -1) +{ + TempLabel = NewLabel(Label_NWP(i,j,k)); + Label_NWP(i,j,k) = TempLabel; + } + } + } + } + //....................................................................... + // Sort WP components by volume + //....................................................................... + OldLabel.resize(NumberComponents_WP); + for (a=0; a -1) +{ + TempLabel = NewLabel(Label_WP(i,j,k)); + Label_WP(i,j,k) = TempLabel; + } + } + } + } + //....................................................................... + +} diff --git a/common/TwoPhase.h b/common/TwoPhase.h index c343c17b..57bd9813 100644 --- a/common/TwoPhase.h +++ b/common/TwoPhase.h @@ -1,4 +1,7 @@ // Header file for two-phase averaging class +#ifndef TwoPhase_INC +#define TwoPhase_INC + #include #include "pmmc.h" @@ -13,42 +16,6 @@ #include "IO/Reader.h" #include "IO/Writer.h" -#define BLOB_AVG_COUNT 33 - -// Array access for averages defined by the following -#define VOL 0 -#define TRIMVOL 1 -#define PRS 2 -#define AWN 3 -#define AWS 4 -#define ANS 5 -#define LWNS 6 -#define JWN 7 -#define KWN 8 -#define CWNS 9 -#define KNWNS 10 -#define KGWNS 11 -#define VX 12 -#define VY 13 -#define VZ 14 -#define VSQ 15 -#define VWNX 16 -#define VWNY 17 -#define VWNZ 18 -#define VWNSX 19 -#define VWNSY 20 -#define VWNSZ 21 -#define GWNXX 22 -#define GWNYY 23 -#define GWNZZ 24 -#define GWNXY 25 -#define GWNXZ 26 -#define GWNYZ 27 -#define TRAWN 28 -#define TRJWN 29 -#define CMX 30 -#define CMY 31 -#define CMZ 32 class TwoPhase{ @@ -144,8 +111,9 @@ public: //........................................................................... int Nx,Ny,Nz; IntArray PhaseID; // Phase ID array (solid=0, non-wetting=1, wetting=2) - IntArray Label_WP; - IntArray Label_NWP; + BlobIDArray Label_WP; // Wetting phase label + BlobIDArray Label_NWP; // Non-wetting phase label index (0:nblobs-1) + std::vector Label_NWP_map; // Non-wetting phase label for each index DoubleArray SDn; DoubleArray SDs; DoubleArray Phase; @@ -169,102 +137,8 @@ public: DoubleArray ComponentAverages_WP; DoubleArray ComponentAverages_NWP; //........................................................................... - TwoPhase(Domain &dm) : Dm(dm){ - Nx=dm.Nx; Ny=dm.Ny; Nz=dm.Nz; - Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz*1.0; - - // Global arrays - PhaseID.resize(Nx,Ny,Nz); - Label_WP.resize(Nx,Ny,Nz); - Label_NWP.resize(Nx,Ny,Nz); - SDn.resize(Nx,Ny,Nz); - SDs.resize(Nx,Ny,Nz); - Phase.resize(Nx,Ny,Nz); - Press.resize(Nx,Ny,Nz); - dPdt.resize(Nx,Ny,Nz); - MeanCurvature.resize(Nx,Ny,Nz); - GaussCurvature.resize(Nx,Ny,Nz); - SDs_x.resize(Nx,Ny,Nz); // Gradient of the signed distance - SDs_y.resize(Nx,Ny,Nz); - SDs_z.resize(Nx,Ny,Nz); - SDn_x.resize(Nx,Ny,Nz); // Gradient of the signed distance - SDn_y.resize(Nx,Ny,Nz); - SDn_z.resize(Nx,Ny,Nz); - DelPhi.resize(Nx,Ny,Nz); - Phase_tplus.resize(Nx,Ny,Nz); - Phase_tminus.resize(Nx,Ny,Nz); - Vel_x.resize(Nx,Ny,Nz); // Gradient of the phase indicator field - Vel_y.resize(Nx,Ny,Nz); - Vel_z.resize(Nx,Ny,Nz); - //......................................... - // Allocate cube storage space - CubeValues.resize(2,2,2); - nw_tris.resize(3,20); - ns_tris.resize(3,20); - ws_tris.resize(3,20); - nws_seg.resize(2,20); - local_sol_tris.resize(3,18); - nw_pts=DTMutableList(20); - ns_pts=DTMutableList(20); - ws_pts=DTMutableList(20); - nws_pts=DTMutableList(20); - local_nws_pts=DTMutableList(20); - local_sol_pts=DTMutableList(20); - tmp=DTMutableList(20); - //......................................... - Values.resize(20); - DistanceValues.resize(20); - KGwns_values.resize(20); - KNwns_values.resize(20); - InterfaceSpeed.resize(20); - NormalVector.resize(60); - //......................................... - van.resize(3); - vaw.resize(3); - vawn.resize(3); - vawns.resize(3); - Gwn.resize(6); - Gns.resize(6); - Gws.resize(6); - van_global.resize(3); - vaw_global.resize(3); - vawn_global.resize(3); - vawns_global.resize(3); - Gwn_global.resize(6); - Gns_global.resize(6); - Gws_global.resize(6); - //......................................... - if (Dm.rank==0){ - TIMELOG = fopen("timelog.tcat","a+"); - if (fseek(TIMELOG,0,SEEK_SET) == fseek(TIMELOG,0,SEEK_CUR)){ - // If timelog is empty, write a short header to list the averages - //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); - fprintf(TIMELOG,"time dEs "); // Timestep, Change in Surface Energy - fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns cwns KNwns KGwns "); // Scalar averages - fprintf(TIMELOG,"vawx vawy vawz vanx vany vanz "); // Velocity averages - fprintf(TIMELOG,"vawnx vawny vawnz vawnsx vawnsy vawnsz "); - fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors - fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz "); - fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz "); - fprintf(TIMELOG,"trawn trJwn trRwn\n"); // trimmed curvature for wn surface - //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); - } - - NWPLOG = fopen("components.NWP.tcat","a+"); - fprintf(NWPLOG,"time label vol pn awn ans Jwn Kwn lwns cwns "); - fprintf(NWPLOG,"vx vy vz vwnx vwny vwnz vwnsx vwnsy vwnsz vsq "); - fprintf(NWPLOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz Cx Cy Cz trawn trJwn\n"); - - WPLOG = fopen("components.WP.tcat","a+"); - fprintf(WPLOG,"time label vol pw awn ans Jwn Kwn lwns cwns "); - fprintf(WPLOG,"vx vy vz vwnx vwny vwnz vwnsx vwnsy vwnsz vsq "); - fprintf(WPLOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz Cx Cy Cz trawn trJwn\n"); - } - } - ~TwoPhase(){ - - } - + TwoPhase(Domain &dm); + ~TwoPhase(); void Initialize(); // void SetupCubes(Domain &Dm); void UpdateMeshValues(); @@ -283,1124 +157,6 @@ public: void PrintComponents(int timestep); }; -void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, DoubleArray &DistData){ -/* double factor,temp,value; - factor=0.5/Beta; - for (int n=0; n 0.999 ) DistData[n] = 4.0; - else if (value < -0.999 ) DistData[n] = -4.0; - else DistData[n] = factor*log((1.0+value)/(1.0-value)); - if (DistData[n] > 1.0) DistData[n] = 1.0; - if (DistData[n] < -1.0) DistData[n] = -1.0; - } - // Initialize to -1,1 (segmentation) - for (int k=0; k 1.0) DistData(i,j,k) = 1.0; - else if (temp < -1.0) DistData(i,j,k) = -1.0; - else DistData(i,j,k) = temp; - } - } - } - - SSO(DistData,Dm.id,Dm,10); -*/ - for (int k=0; k 0 && Dm.kproc == 0) kmin=4; - if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4; - - for (k=kmin; k 0 ){ - nwp_volume += 0.125; - // velocity - van(0) += 0.125*Vel_x(n); - van(1) += 0.125*Vel_y(n); - van(2) += 0.125*Vel_z(n); - // volume the excludes the interfacial region - if (DelPhi(n) < 1e-4){ - vol_n += 0.125; - // pressure - pan += 0.125*Press(n); - - } - } - else{ - wp_volume += 0.125; - // velocity - vaw(0) += 0.125*Vel_x(n); - vaw(1) += 0.125*Vel_y(n); - vaw(2) += 0.125*Vel_z(n); - if (DelPhi(n) < 1e-4){ - // volume the excludes the interfacial region - vol_w += 0.125; - // pressure - paw += 0.125*Press(n); - - } - } - } - } - - //........................................................................... - // Construct the interfaces and common curve - pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue, - nw_pts, nw_tris, Values, ns_pts, ns_tris, ws_pts, ws_tris, - local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris, - n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris, - n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg, - i, j, k, Nx, Ny, Nz); - - // wn interface averages - if (n_nw_pts > 0){ - awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris); - Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); - Kwn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); - - // Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels) - pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues, - i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn); - - pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues, - i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn); - - // Compute the normal speed of the interface - pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris, - NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris); - } - // wns common curve averages - if (n_local_nws_pts > 0){ - efawns += pmmc_CubeContactAngle(CubeValues,Values,SDn_x,SDn_y,SDn_z,SDs_x,SDs_y,SDs_z, - local_nws_pts,i,j,k,n_local_nws_pts); - - pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z, - local_nws_pts,i,j,k,n_local_nws_pts); - - pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y, - SDs_z, KNwns_values, KGwns_values, KNwns, KGwns, - nws_pts, n_nws_pts, i, j, k); - - lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts); - } - - // Solid interface averagees - if (n_local_sol_tris > 0){ - As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris); - - // Compute the surface orientation and the interfacial area - ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris); - aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris); - } - //........................................................................... - } - } - } -} -void TwoPhase::AssignComponentLabels(){ - - int LabelNWP=1; - //int LabelWP=2; - // NOTE: labeling the wetting phase components is tricky! One sandstone media had over 800,000 components - //NumberComponents_WP = ComputeGlobalPhaseComponent(Dm.Nx-2,Dm.Ny-2,Dm.Nz-2,Dm.rank_info,PhaseID,LabelWP,Label_WP); - // treat all wetting phase is connected - NumberComponents_WP=1; - for (int k=0; k 0 && Dm.kproc == 0) kmin=4; - if (Dm.BoundaryCondition > 0 && Dm.kproc == Dm.nprocz-1) kmax=Nz-4; - - for (k=kmin; k 0.0 && !(LabelNWP < 0) ){ - // volume - ComponentAverages_NWP(VOL,LabelNWP) += 0.125; - // velocity - ComponentAverages_NWP(VX,LabelNWP) += 0.125*Vel_x(n); - ComponentAverages_NWP(VY,LabelNWP) += 0.125*Vel_y(n); - ComponentAverages_NWP(VZ,LabelNWP) += 0.125*Vel_z(n); - // center of mass - ComponentAverages_NWP(CMX,LabelNWP) += 0.125*(i+cube[p][0]+Dm.iproc*Nx); - ComponentAverages_NWP(CMY,LabelNWP) += 0.125*(j+cube[p][1]+Dm.jproc*Ny); - ComponentAverages_NWP(CMZ,LabelNWP) += 0.125*(k+cube[p][2]+Dm.kproc*Nz); - - // twice the kinetic energy - ComponentAverages_NWP(VSQ,LabelNWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n)); - - // volume the for pressure averaging excludes the interfacial region - if (DelPhi(n) < 1e-4 ){ - ComponentAverages_NWP(TRIMVOL,LabelNWP) += 0.125; - ComponentAverages_NWP(PRS,LabelNWP) += 0.125*Press(n); - } - } - else if (!(LabelWP < 0)){ - ComponentAverages_WP(VOL,LabelWP) += 0.125; - // velocity - ComponentAverages_WP(VX,LabelWP) += 0.125*Vel_x(n); - ComponentAverages_WP(VY,LabelWP)+= 0.125*Vel_y(n); - ComponentAverages_WP(VZ,LabelWP) += 0.125*Vel_z(n); - // Center of mass - ComponentAverages_WP(CMX,LabelWP) += 0.125*(i+cube[p][0]+Dm.iproc*Nx); - ComponentAverages_WP(CMY,LabelWP) += 0.125*(j+cube[p][1]+Dm.jproc*Ny); - ComponentAverages_WP(CMZ,LabelWP) += 0.125*(k+cube[p][2]+Dm.kproc*Nz); - // twice the kinetic energy - ComponentAverages_WP(VSQ,LabelWP) += 0.125*(Vel_x(n)*Vel_x(n)+Vel_y(n)*Vel_y(n)+Vel_z(n)*Vel_z(n)); - - // volume the for pressure averaging excludes the interfacial region - if (DelPhi(n) < 1e-4){ - ComponentAverages_WP(TRIMVOL,LabelWP) += 0.125; - ComponentAverages_WP(PRS,LabelWP) += 0.125*Press(n); - } - } - } - } - //........................................................................... - // Construct the interfaces and common curve - pmmc_ConstructLocalCube(SDs, SDn, solid_isovalue, fluid_isovalue, - nw_pts, nw_tris, Values, ns_pts, ns_tris, ws_pts, ws_tris, - local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris, - n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris, - n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg, - i, j, k, Nx, Ny, Nz); - - //........................................................................... - // wn interface averages - if (n_nw_pts>0 && LabelNWP >=0 && LabelWP >=0 ){ - // Mean curvature - TempLocal = pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); - ComponentAverages_WP(JWN,LabelWP) += TempLocal; - ComponentAverages_NWP(JWN,LabelNWP) += TempLocal; - - // Trimmed Mean curvature - pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues, - i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn); - ComponentAverages_WP(TRAWN,LabelWP) += trawn; - ComponentAverages_WP(TRJWN,LabelWP) += trJwn; - ComponentAverages_NWP(TRAWN,LabelNWP) += trawn; - ComponentAverages_NWP(TRJWN,LabelNWP) += trJwn; - - // Gaussian curvature - TempLocal = pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); - ComponentAverages_WP(KWN,LabelWP) += TempLocal; - ComponentAverages_NWP(KWN,LabelNWP) += TempLocal; - - // Compute the normal speed of the interface - pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris, - NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris); - ComponentAverages_WP(VWNX,LabelWP) += vawn(0); - ComponentAverages_WP(VWNY,LabelWP) += vawn(1); - ComponentAverages_WP(VWNZ,LabelWP) += vawn(2); - ComponentAverages_NWP(VWNX,LabelNWP) += vawn(0); - ComponentAverages_NWP(VWNY,LabelNWP) += vawn(1); - ComponentAverages_NWP(VWNZ,LabelNWP) += vawn(2); - - // Interfacial Area - TempLocal = pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris); - ComponentAverages_WP(AWN,LabelWP) += TempLocal; - ComponentAverages_NWP(AWN,LabelNWP) += TempLocal; - - ComponentAverages_WP(GWNXX,LabelWP) += Gwn(0); - ComponentAverages_WP(GWNYY,LabelWP) += Gwn(1); - ComponentAverages_WP(GWNZZ,LabelWP) += Gwn(2); - ComponentAverages_WP(GWNXY,LabelWP) += Gwn(3); - ComponentAverages_WP(GWNXZ,LabelWP) += Gwn(4); - ComponentAverages_WP(GWNYZ,LabelWP) += Gwn(5); - - ComponentAverages_NWP(GWNXX,LabelNWP) += Gwn(0); - ComponentAverages_NWP(GWNYY,LabelNWP) += Gwn(1); - ComponentAverages_NWP(GWNZZ,LabelNWP) += Gwn(2); - ComponentAverages_NWP(GWNXY,LabelNWP) += Gwn(3); - ComponentAverages_NWP(GWNXZ,LabelNWP) += Gwn(4); - ComponentAverages_NWP(GWNYZ,LabelNWP) += Gwn(5); - - } - //........................................................................... - // Common curve averages - if (n_local_nws_pts > 0 && LabelNWP >=0 && LabelWP >=0){ - // Contact angle - TempLocal = pmmc_CubeContactAngle(CubeValues,Values,SDn_x,SDn_y,SDn_z,SDs_x,SDs_y,SDs_z, - local_nws_pts,i,j,k,n_local_nws_pts); - ComponentAverages_WP(CWNS,LabelWP) += TempLocal; - ComponentAverages_NWP(CWNS,LabelNWP) += TempLocal; - - // Kinematic velocity of the common curve - pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z, - local_nws_pts,i,j,k,n_local_nws_pts); - ComponentAverages_WP(VWNSX,LabelWP) += vawns(0); - ComponentAverages_WP(VWNSY,LabelWP) += vawns(1); - ComponentAverages_WP(VWNSZ,LabelWP) += vawns(2); - ComponentAverages_NWP(VWNSX,LabelNWP) += vawns(0); - ComponentAverages_NWP(VWNSY,LabelNWP) += vawns(1); - ComponentAverages_NWP(VWNSZ,LabelNWP) += vawns(2); - - // Curvature of the common curve - pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y, - SDs_z, KNwns_values, KGwns_values, KNwns, KGwns, - nws_pts, n_nws_pts, i, j, k); - ComponentAverages_WP(KNWNS,LabelWP) += KNwns; - ComponentAverages_WP(KGWNS,LabelWP) += KGwns; - ComponentAverages_NWP(KNWNS,LabelNWP) += KNwns; - ComponentAverages_NWP(KGWNS,LabelNWP) += KGwns; - - // Length of the common curve - TempLocal = pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts); - ComponentAverages_NWP(LWNS,LabelNWP) += TempLocal; - } - //........................................................................... - // Solid interface averages - if (n_local_sol_pts > 0 && LabelWP >=0){ - As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris); - // Compute the surface orientation and the interfacial area - - TempLocal = pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris); - ComponentAverages_WP(AWS,LabelWP) += TempLocal; - } - if (n_ns_pts > 0 && LabelNWP >=0 ){ - TempLocal = pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris); - ComponentAverages_NWP(ANS,LabelNWP) += TempLocal; - } - //........................................................................... - - } - } - } - - if (Dm.rank==0){ - printf("Component averages computed locally -- reducing result... \n"); - } - // Globally reduce the non-wetting phase averages - RecvBuffer.resize(BLOB_AVG_COUNT,NumberComponents_NWP); - -/* for (int b=0; b 0.0){ - double Vn,pn,awn,ans,Jwn,Kwn,lwns,cwns,vsq; - - Vn = ComponentAverages_NWP(VOL,b); - awn = ComponentAverages_NWP(AWN,b); - ans = ComponentAverages_NWP(ANS,b); - van(0) = ComponentAverages_NWP(VX,b)/Vn; - van(1) = ComponentAverages_NWP(VY,b)/Vn; - van(2) = ComponentAverages_NWP(VZ,b)/Vn; - vsq = ComponentAverages_NWP(VSQ,b)/Vn; - NULL_USE(ans); - - if (ComponentAverages_NWP(TRIMVOL,b) > 0.0){ - pn = ComponentAverages_NWP(PRS,b)/ComponentAverages_NWP(TRIMVOL,b); - } - else pn = 0.0; - - if (awn != 0.0){ - Jwn = ComponentAverages_NWP(JWN,b)/awn; - Kwn = ComponentAverages_NWP(KWN,b)/awn; - vawn(0) = ComponentAverages_NWP(VWNSX,b)/awn; - vawn(1) = ComponentAverages_NWP(VWNSY,b)/awn; - vawn(2) = ComponentAverages_NWP(VWNSZ,b)/awn; - Gwn(0) = ComponentAverages_NWP(GWNXX,b)/awn; - Gwn(1) = ComponentAverages_NWP(GWNYY,b)/awn; - Gwn(2) = ComponentAverages_NWP(GWNZZ,b)/awn; - Gwn(3) = ComponentAverages_NWP(GWNXY,b)/awn; - Gwn(4) = ComponentAverages_NWP(GWNXZ,b)/awn; - Gwn(5) = ComponentAverages_NWP(GWNYZ,b)/awn; - } - else Jwn=Kwn=0.0; - - trawn = ComponentAverages_NWP(TRAWN,b); - trJwn = ComponentAverages_NWP(TRJWN,b); - if (trawn > 0.0) trJwn /= trawn; - - lwns = ComponentAverages_NWP(LWNS,b); - if (lwns != 0.0){ - cwns = ComponentAverages_NWP(CWNS,b)/lwns; - vawns(0) = ComponentAverages_NWP(VWNSX,b)/lwns; - vawns(1) = ComponentAverages_NWP(VWNSY,b)/lwns; - vawns(2) = ComponentAverages_NWP(VWNSZ,b)/lwns; - } - else cwns=0.0; - - ComponentAverages_NWP(PRS,b) = pn; - ComponentAverages_NWP(VX,b) = van(0); - ComponentAverages_NWP(VY,b) = van(1); - ComponentAverages_NWP(VZ,b) = van(2); - ComponentAverages_NWP(VSQ,b) = vsq; - - ComponentAverages_NWP(JWN,b) = Jwn; - ComponentAverages_NWP(KWN,b) = Kwn; - ComponentAverages_NWP(VWNX,b) = vawn(0); - ComponentAverages_NWP(VWNY,b) = vawn(1); - ComponentAverages_NWP(VWNZ,b) = vawn(2); - - ComponentAverages_NWP(GWNXX,b) = Gwn(0); - ComponentAverages_NWP(GWNYY,b) = Gwn(1); - ComponentAverages_NWP(GWNZZ,b) = Gwn(2); - ComponentAverages_NWP(GWNXY,b) = Gwn(3); - ComponentAverages_NWP(GWNXZ,b) = Gwn(4); - ComponentAverages_NWP(GWNYZ,b) = Gwn(5); - - ComponentAverages_NWP(CWNS,b) = cwns; - ComponentAverages_NWP(VWNSX,b) = vawns(0); - ComponentAverages_NWP(VWNSY,b) = vawns(1); - ComponentAverages_NWP(VWNSZ,b) = vawns(2); - - ComponentAverages_NWP(CMX,b) /= Vn; - ComponentAverages_NWP(CMY,b) /= Vn; - ComponentAverages_NWP(CMZ,b) /= Vn; - - ComponentAverages_NWP(TRJWN,b) = trJwn; - - } - } - - // reduce the wetting phase averages - for (int b=0; b 0.0){ - double Vw,pw,awn,ans,Jwn,Kwn,lwns,cwns,vsq; - Vw = ComponentAverages_WP(VOL,b); - awn = ComponentAverages_WP(AWN,b); - ans = ComponentAverages_WP(ANS,b); - vaw(0) = ComponentAverages_WP(VX,b)/Vw; - vaw(1) = ComponentAverages_WP(VY,b)/Vw; - vaw(2) = ComponentAverages_WP(VZ,b)/Vw; - vsq = ComponentAverages_WP(VSQ,b)/Vw; - NULL_USE(ans); - - if (ComponentAverages_WP(TRIMVOL,b) > 0.0){ - pw = ComponentAverages_WP(PRS,b)/ComponentAverages_WP(TRIMVOL,b); - } - else pw = 0.0; - - if (awn != 0.0){ - Jwn = ComponentAverages_WP(JWN,b)/awn; - Kwn = ComponentAverages_WP(KWN,b)/awn; - vawn(0) = ComponentAverages_WP(VWNSX,b)/awn; - vawn(1) = ComponentAverages_WP(VWNSY,b)/awn; - vawn(2) = ComponentAverages_WP(VWNSZ,b)/awn; - Gwn(0) = ComponentAverages_WP(GWNXX,b)/awn; - Gwn(1) = ComponentAverages_WP(GWNYY,b)/awn; - Gwn(2) = ComponentAverages_WP(GWNZZ,b)/awn; - Gwn(3) = ComponentAverages_WP(GWNXY,b)/awn; - Gwn(4) = ComponentAverages_WP(GWNXZ,b)/awn; - Gwn(5) = ComponentAverages_WP(GWNYZ,b)/awn; - } - else Jwn=Kwn=0.0; - - trawn = ComponentAverages_WP(TRAWN,b); - trJwn = ComponentAverages_WP(TRJWN,b); - if (trawn > 0.0) trJwn /= trawn; - - lwns = ComponentAverages_WP(LWNS,b); - if (lwns != 0.0){ - cwns = ComponentAverages_WP(CWNS,b)/lwns; - vawns(0) = ComponentAverages_WP(VWNSX,b)/lwns; - vawns(1) = ComponentAverages_WP(VWNSY,b)/lwns; - vawns(2) = ComponentAverages_WP(VWNSZ,b)/lwns; - } - else cwns=0.0; - - ComponentAverages_WP(PRS,b) = pw; - ComponentAverages_WP(VX,b) = vaw(0); - ComponentAverages_WP(VY,b) = vaw(1); - ComponentAverages_WP(VZ,b) = vaw(2); - ComponentAverages_WP(VSQ,b) = vsq; - - ComponentAverages_WP(JWN,b) = Jwn; - ComponentAverages_WP(KWN,b) = Kwn; - ComponentAverages_WP(VWNX,b) = vawn(0); - ComponentAverages_WP(VWNY,b) = vawn(1); - ComponentAverages_WP(VWNZ,b) = vawn(2); - - ComponentAverages_WP(GWNXX,b) = Gwn(0); - ComponentAverages_WP(GWNYY,b) = Gwn(1); - ComponentAverages_WP(GWNZZ,b) = Gwn(2); - ComponentAverages_WP(GWNXY,b) = Gwn(3); - ComponentAverages_WP(GWNXZ,b) = Gwn(4); - ComponentAverages_WP(GWNYZ,b) = Gwn(5); - - ComponentAverages_WP(CWNS,b) = cwns; - ComponentAverages_WP(VWNSX,b) = vawns(0); - ComponentAverages_WP(VWNSY,b) = vawns(1); - ComponentAverages_WP(VWNSZ,b) = vawns(2); - - ComponentAverages_WP(TRJWN,b) = trJwn; - } - } -} - -void TwoPhase::WriteSurfaces(int logcount){ - - int i,j,k; - int ncubes=(Nx-1)*(Ny-1)*(Nz-1); - Point P,A,B,C; - - std::shared_ptr wn_mesh( new IO::TriList() ); - wn_mesh->A.reserve(8*ncubes); - wn_mesh->B.reserve(8*ncubes); - wn_mesh->C.reserve(8*ncubes); - - std::shared_ptr ns_mesh( new IO::TriList() ); - ns_mesh->A.reserve(8*ncubes); - ns_mesh->B.reserve(8*ncubes); - ns_mesh->C.reserve(8*ncubes); - - std::shared_ptr ws_mesh( new IO::TriList() ); - ws_mesh->A.reserve(8*ncubes); - ws_mesh->B.reserve(8*ncubes); - ws_mesh->C.reserve(8*ncubes); - - for (k=1; kA.push_back(A); - wn_mesh->B.push_back(B); - wn_mesh->C.push_back(C); - } - for (int r=0;rA.push_back(A); - ws_mesh->B.push_back(B); - ws_mesh->C.push_back(C); - } - for (int r=0;rA.push_back(A); - ns_mesh->B.push_back(B); - ns_mesh->C.push_back(C); - } - } - } - } - - std::vector meshData(4); - meshData[0].meshName = "wn-tris"; - meshData[0].mesh = wn_mesh; - meshData[1].meshName = "ws-tris"; - meshData[1].mesh = ws_mesh; - meshData[2].meshName = "ns-tris"; - meshData[2].mesh = ns_mesh; - IO::writeData( logcount, meshData, 2); - -} - -void TwoPhase::Reduce(){ - int i; - double iVol_global=1.0/Volume; - //........................................................................... - MPI_Barrier(Dm.Comm); - MPI_Allreduce(&nwp_volume,&nwp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&wp_volume,&wp_volume_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&awn,&awn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&ans,&ans_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&aws,&aws_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&lwns,&lwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&As,&As_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&Jwn,&Jwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&Kwn,&Kwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&KGwns,&KGwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&KNwns,&KNwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - // Phase averages - MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&paw,&paw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&pan,&pan_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&vaw(0),&vaw_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&van(0),&van_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&vawn(0),&vawn_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&vawns(0),&vawns_global(0),3,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&Gwn(0),&Gwn_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&Gns(0),&Gns_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&Gws(0),&Gws_global(0),6,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&trawn,&trawn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&trJwn,&trJwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Allreduce(&trRwn,&trRwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); - MPI_Barrier(Dm.Comm); - - // Normalize the phase averages - // (density of both components = 1.0) - if (vol_w_global > 0.0){ - paw_global = paw_global / vol_w_global; - } - if (wp_volume_global > 0.0){ - vaw_global(0) = vaw_global(0) / wp_volume_global; - vaw_global(1) = vaw_global(1) / wp_volume_global; - vaw_global(2) = vaw_global(2) / wp_volume_global; - - } - if (vol_n_global > 0.0){ - pan_global = pan_global / vol_n_global; - } - if (nwp_volume_global > 0.0){ - van_global(0) = van_global(0) / nwp_volume_global; - van_global(1) = van_global(1) / nwp_volume_global; - van_global(2) = van_global(2) / nwp_volume_global; - } - // Normalize surface averages by the interfacial area - if (awn_global > 0.0){ - Jwn_global /= awn_global; - Kwn_global /= awn_global; - for (i=0; i<3; i++) vawn_global(i) /= awn_global; - for (i=0; i<6; i++) Gwn_global(i) /= awn_global; - } - if (lwns_global > 0.0){ - efawns_global /= lwns_global; - KNwns_global /= lwns_global; - KGwns_global /= lwns_global; - for (i=0; i<3; i++) vawns_global(i) /= lwns_global; - } - if (trawn_global > 0.0){ - trJwn_global /= trawn_global; - trRwn_global /= trawn_global; - trRwn_global = 2.0*fabs(trRwn_global); - trJwn_global = fabs(trJwn_global); - } - - if (ans_global > 0.0) for (i=0; i<6; i++) Gns_global(i) /= ans_global; - if (aws_global > 0.0) for (i=0; i<6; i++) Gws_global(i) /= aws_global; - - //sat_w = 1.0 - nwp_volume_global*iVol_global/porosity; - sat_w = 1.0 - nwp_volume_global/(nwp_volume_global+wp_volume_global); - // Compute the specific interfacial areas and common line length (dimensionless per unit volume) - awn_global = awn_global*iVol_global; - ans_global = ans_global*iVol_global; - aws_global = aws_global*iVol_global; - dEs = dEs*iVol_global; - lwns_global = lwns_global*iVol_global; -} - -void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT){ - awn_global *= D; - ans_global *= D; - ans_global *= D; - lwns_global *= D*D; -} - -void TwoPhase::PrintAll(int timestep){ - if (Dm.rank==0){ - fprintf(TIMELOG,"%i %.5g ",timestep-5,dEs); // change in surface energy - fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure - fprintf(TIMELOG,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas - fprintf(TIMELOG,"%.5g %.5g ",Jwn_global, Kwn_global); // curvature of wn interface - fprintf(TIMELOG,"%.5g ",lwns_global); // common curve length - fprintf(TIMELOG,"%.5g ",efawns_global); // average contact angle - fprintf(TIMELOG,"%.5g %.5g ",KNwns_global, KGwns_global); // curvature of wn interface - fprintf(TIMELOG,"%.5g %.5g %.5g ",vaw_global(0),vaw_global(1),vaw_global(2)); // average velocity of w phase - fprintf(TIMELOG,"%.5g %.5g %.5g ",van_global(0),van_global(1),van_global(2)); // average velocity of n phase - fprintf(TIMELOG,"%.5g %.5g %.5g ",vawn_global(0),vawn_global(1),vawn_global(2)); // velocity of wn interface - fprintf(TIMELOG,"%.5g %.5g %.5g ",vawns_global(0),vawns_global(1),vawns_global(2)); // velocity of wn interface - fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", - Gwn_global(0),Gwn_global(1),Gwn_global(2),Gwn_global(3),Gwn_global(4),Gwn_global(5)); // orientation of wn interface - fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", - Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5)); // orientation of ns interface - fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", - Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface - fprintf(TIMELOG,"%.5g %.5g %.5g\n",trawn_global, trJwn_global, trRwn_global); // Trimmed curvature - fflush(TIMELOG); - } -} - -void TwoPhase::PrintComponents(int timestep){ - if (Dm.rank==0){ - printf("PRINT %i COMPONENT AVEREAGES: time = %i \n",(int)ComponentAverages_NWP.size(1),timestep); - for (int b=0; 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 ",ComponentAverages_NWP(GWNYZ,b)); - fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CMX,b)); - fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CMY,b)); - fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(CMZ,b)); - fprintf(NWPLOG,"%.5g ",ComponentAverages_NWP(TRAWN,b)); - fprintf(NWPLOG,"%.5g\n",ComponentAverages_NWP(TRJWN,b)); - } - } - fflush(NWPLOG); - - for (int b=0; 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 ",ComponentAverages_WP(GWNYZ,b)); - fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CMX,b)); - fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CMY,b)); - fprintf(WPLOG,"%.5g ",ComponentAverages_WP(CMZ,b)); - fprintf(WPLOG,"%.5g ",ComponentAverages_WP(TRAWN,b)); - fprintf(WPLOG,"%.5g\n",ComponentAverages_WP(TRJWN,b)); - - } - } - fflush(WPLOG); - } -} - -inline int TwoPhase::GetCubeLabel(int i, int j, int k, IntArray &BlobLabel){ - int label; - label=BlobLabel(i,j,k); - label=max(label,BlobLabel(i+1,j,k)); - label=max(label,BlobLabel(i,j+1,k)); - label=max(label,BlobLabel(i+1,j+1,k)); - label=max(label,BlobLabel(i,j,k+1)); - label=max(label,BlobLabel(i+1,j,k+1)); - label=max(label,BlobLabel(i,j+1,k+1)); - label=max(label,BlobLabel(i+1,j+1,k+1)); - - return label; -} - -void TwoPhase::SortBlobs(){ - //printf("Sorting the blobs based on volume \n"); - //printf("-----------------------------------------------\n"); - int TempLabel,a,aa,bb,i,j,k,idx; - double TempValue; - //....................................................................... - // Sort NWP components by volume - //....................................................................... - IntArray OldLabel(NumberComponents_NWP); - for (a=0; a -1){ - TempLabel = NewLabel(Label_NWP(i,j,k)); - Label_NWP(i,j,k) = TempLabel; - } - } - } - } - //....................................................................... - // Sort WP components by volume - //....................................................................... - OldLabel.resize(NumberComponents_WP); - for (a=0; a -1){ - TempLabel = NewLabel(Label_WP(i,j,k)); - Label_WP(i,j,k) = TempLabel; - } - } - } - } - //....................................................................... - -} diff --git a/common/pmmc.h b/common/pmmc.h index 8aefd466..b42725cd 100644 --- a/common/pmmc.h +++ b/common/pmmc.h @@ -1,3 +1,7 @@ +#ifndef pmmc_INC +#define pmmc_INC + + #include #include #include @@ -10,7 +14,7 @@ using namespace std; -int edgeTable[256]={ +static int edgeTable[256]={ 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, @@ -43,8 +47,8 @@ int edgeTable[256]={ 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 }; -char triTable[256][16] = -{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, +static char triTable[256][16] = + {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, @@ -735,7 +739,7 @@ inline bool Solid( DoubleArray &A, int i, int j, int k){ return X; } //------------------------------------------------------------------------------- -Point VertexInterp(const Point &p1, const Point &p2, double valp1, double valp2) +inline Point VertexInterp(const Point &p1, const Point &p2, double valp1, double valp2) { return (p1 + (-valp1 / (valp2 - valp1)) * (p2 - p1)); } @@ -4032,3 +4036,7 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray //............................................................................. } //-------------------------------------------------------------------------------------------------------- + + +#endif + diff --git a/pmmc/pmmc.cpp b/pmmc/pmmc.cpp index 9a332598..5fd4d745 100644 --- a/pmmc/pmmc.cpp +++ b/pmmc/pmmc.cpp @@ -11,7 +11,7 @@ using namespace std; //-------------------------------------------------------------------------------------------------------- -inline int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator, +int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indicator, DoubleArray &F, DoubleArray &S, double vf, double vs, int startx, int starty, int startz, IntArray &temp) { @@ -140,7 +140,7 @@ inline int ComputeBlob(IntArray &blobs, int &nblobs, int &ncubes, IntArray &indi } //-------------------------------------------------------------------------------------------------- /* -inline DoubleArray SOLVE( DoubleArray &A, DoubleArray &b) +DoubleArray SOLVE( DoubleArray &A, DoubleArray &b) { // solves the system A*x = b exactly @@ -158,7 +158,7 @@ inline DoubleArray SOLVE( DoubleArray &A, DoubleArray &b) return solution; } //----------------------------------------------------------------------------- -inline Point NEWTON(Point x0, DoubleArray &F, double &v,DoubleArray &S, int i, int j, int k, int &newton_steps) +Point NEWTON(Point x0, DoubleArray &F, double &v,DoubleArray &S, int i, int j, int k, int &newton_steps) { Point pt; // Performs a Newton iteration to compute the point x from initial guess x0 @@ -321,7 +321,7 @@ inline Point NEWTON(Point x0, DoubleArray &F, double &v,DoubleArray &S, int i, */ //-------------------------------------------------------------------------------------------------- -inline bool vertexcheck(Point &P, int n, int pos, DTMutableList &cellvertices){ +bool vertexcheck(Point &P, int n, int pos, DTMutableList &cellvertices){ // returns true if P is a new vertex (one previously unencountered bool V = 1; @@ -337,7 +337,7 @@ inline bool vertexcheck(Point &P, int n, int pos, DTMutableList &cellvert //-------------------------------------------------------------------------------------------------- -inline bool ShareSide( Point &A, Point &B) +bool ShareSide( Point &A, Point &B) { // returns true if points A and B share an x,y, or z coordinate bool l = 0; @@ -360,7 +360,7 @@ inline bool ShareSide( Point &A, Point &B) } //-------------------------------------------------------------------------------------------------- -inline bool Interface( DoubleArray &A, const double v, int i, int j, int k){ +bool Interface( DoubleArray &A, const double v, int i, int j, int k){ // returns true if grid cell i, j, k contains a section of the interface bool Y = 0; @@ -416,7 +416,7 @@ inline bool Interface( DoubleArray &A, const double v, int i, int j, int k){ } //-------------------------------------------------------------------------------------------------- -inline bool Fluid_Interface( DoubleArray &A, DoubleArray &S, const double v, int i, int j, int k){ +bool Fluid_Interface( DoubleArray &A, DoubleArray &S, const double v, int i, int j, int k){ // returns true if grid cell i, j, k contains a section of the interface bool Y = 0; @@ -470,7 +470,7 @@ inline bool Fluid_Interface( DoubleArray &A, DoubleArray &S, const double v, in return Y; } //-------------------------------------------------------------------------------------------------- -inline bool Solid( DoubleArray &A, int i, int j, int k){ +bool Solid( DoubleArray &A, int i, int j, int k){ bool X = 0; @@ -504,7 +504,7 @@ inline bool Solid( DoubleArray &A, int i, int j, int k){ } //------------------------------------------------------------------------------- //------------------------------------------------------------------------------- -inline void SOL_SURF(DoubleArray &A, const double &v, DoubleArray &B, const double &isovalue, +void SOL_SURF(DoubleArray &A, const double &v, DoubleArray &B, const double &isovalue, int i,int j,int k, int m, int n, int o, DTMutableList &cellvertices, int &lengthvertices, IntArray &Tlist, int &nTris, DoubleArray &values){ @@ -941,7 +941,7 @@ inline void SOL_SURF(DoubleArray &A, const double &v, DoubleArray &B, const doub } } //------------------------------------------------------------------------------- -inline void TRIM(DTMutableList &local_sol_pts, int &n_local_sol_pts, double isovalue, +void TRIM(DTMutableList &local_sol_pts, int &n_local_sol_pts, double isovalue, IntArray &local_sol_tris, int &n_local_sol_tris, DTMutableList &ns_pts, int &n_ns_pts, IntArray &ns_tris, int &n_ns_tris, DTMutableList &ws_pts, int &n_ws_pts, @@ -1374,7 +1374,7 @@ inline void TRIM(DTMutableList &local_sol_pts, int &n_local_sol_pts, doub } } //------------------------------------------------------------------------------- -inline void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k, +void MC( DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k, DTMutableList &nw_pts, int &n_nw_pts, IntArray &nw_tris, int &n_nw_tris) { @@ -2128,7 +2128,7 @@ else if ( A(i,j+1,k+1) == 0){ } //------------------------------------------------------------------------------- -inline void EDGE(DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k, int &m, int &n, int &o, +void EDGE(DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, int &k, int &m, int &n, int &o, DTMutableList &nw_pts, int &n_nw_pts, IntArray &nw_tris, int &n_nw_tris, DTMutableList &local_nws_pts, int &n_local_nws_pts) { @@ -2491,7 +2491,7 @@ inline void EDGE(DoubleArray &A, double &v, DoubleArray &solid, int &i, int &j, } //-------------------------------------------------------------------------------------------------------- -inline void ComputeAreasPMMC(IntArray &cubeList, int start, int finish, +void ComputeAreasPMMC(IntArray &cubeList, int start, int finish, DoubleArray &F, DoubleArray &S, double vF, double vS, double &blob_volume, double &ans, double &aws, double &awn, double &lwns, int Nx, int Ny, int Nz) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index fe187fa5..d349fb8c 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -27,9 +27,9 @@ ADD_LBPM_TEST( TestInterfaceSpeed ) ADD_LBPM_TEST( TestSphereCurvature ) ADD_LBPM_TEST( TestContactAngle ) ADD_LBPM_TEST_1_2_4( TestTwoPhase ) +ADD_LBPM_TEST_1_2_4( TestBlobIdentify ) ADD_LBPM_TEST_PARALLEL( TestTwoPhase 8 ) ADD_LBPM_TEST_PARALLEL( TestBlobAnalyze 8 ) -ADD_LBPM_TEST_PARALLEL( TestBlobIdentify 1 ) ADD_LBPM_TEST_PARALLEL( TestSegDist 8 ) ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 ) ADD_LBPM_TEST_PARALLEL( TestMassConservationD3Q7 1 ) diff --git a/tests/TestBlobIdentify.cpp b/tests/TestBlobIdentify.cpp index 190e1b94..91e75ce2 100644 --- a/tests/TestBlobIdentify.cpp +++ b/tests/TestBlobIdentify.cpp @@ -22,6 +22,16 @@ inline double rand2() } +// Test if all ranks agree on a value +bool allAgree( int x ) { + int min, max; + MPI_Allreduce(&x,&min,1,MPI_INT,MPI_MIN,MPI_COMM_WORLD); + MPI_Allreduce(&x,&max,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD); + return min==max; +} + + +// Structure to hold a bubble struct bubble_struct { Point center; double radius; @@ -209,6 +219,7 @@ int main(int argc, char **argv) fillData.copy(SignDist,SignDistVar->data); fillData.copy(GlobalBlobID,BlobIDVar->data); IO::writeData( 0, meshData, 2 ); + writeIDMap(ID_map_struct(nblobs),0,"lbpm_id_map.txt"); int save_it = 1; // Check the results @@ -222,6 +233,7 @@ int main(int argc, char **argv) // 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"); } + int id_max = nblobs-1; for (int i=0; i<20; i++, save_it++) { // Shift all the data shift_data( Phase, 3, -2, 1, rank_info ); @@ -236,6 +248,11 @@ int main(int argc, char **argv) } // Identify the blob maps and renumber the ids ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2); + std::swap(GlobalBlobID,GlobalBlobID2); + std::vector new_list; + getNewIDs(map,id_max,new_list); + renumberIDs(new_list,GlobalBlobID); + writeIDMap(map,save_it,"lbpm_id_map.txt"); bool pass = (int)map.src_dst.size()==nblobs; pass = pass && map.created.empty(); pass = pass && map.destroyed.empty(); @@ -244,11 +261,11 @@ int main(int argc, char **argv) pass = pass && map.split.empty(); pass = pass && map.merge.empty(); pass = pass && map.merge_split.empty(); + pass = pass && id_max==nblobs-1; 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); @@ -279,6 +296,7 @@ int main(int argc, char **argv) fillData.copy(GlobalBlobID,BlobIDVar->data); IO::writeData( save_it, meshData, 2 ); save_it++; + id_max = nblobs-1; for (int i=0; i<25; i++, save_it++) { // Move the bubbles for (size_t j=0; j new_list; + getNewIDs(map,id_max,new_list); + renumberIDs(new_list,GlobalBlobID); + writeIDMap(map,save_it,"lbpm_id_map.txt"); + if ( !allAgree(id_max) ) { + if ( rank==0 ) + printf("All ranks do not agree on id_max\n"); + N_errors++; + break; + } int N1 = 0; int N2 = 0; if ( rank==0 ) { @@ -350,7 +379,6 @@ int main(int argc, char **argv) N_errors++; } nblobs = nblobs2; - GlobalBlobID = GlobalBlobID2; // Save the results fillData.copy(Phase,PhaseVar->data); fillData.copy(SignDist,SignDistVar->data); diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index 1bf4c66a..21ca9d32 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -677,8 +677,10 @@ int main(int argc, char **argv) ThreadPool::setProcessAffinity(procs); int timestep = -1; AnalysisWaitIdStruct work_ids; - ThreadPool tpool(2); - BlobIDstruct last_ids; + ThreadPool tpool(0); + BlobIDstruct last_ids, last_index; + BlobIDList last_id_map; + writeIDMap(ID_map_struct(),0,id_map_filename); while (timestep < timestepMax && err > tol ) { PROFILE_START("Update"); @@ -784,7 +786,7 @@ int main(int argc, char **argv) timestep++; // Run the analysis, blob identification, and write restart files - run_analysis(timestep,RESTART_INTERVAL,rank_info,Averages,last_ids, + run_analysis(timestep,RESTART_INTERVAL,rank_info,Averages,last_ids,last_index,last_id_map, Nx,Ny,Nz,pBC,beta,err,Phi,Pressure,Velocity,ID,f_even,f_odd,Den, LocalRestartFile,tpool,work_ids); diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h index 769a38ec..baaca18e 100644 --- a/tests/lbpm_color_simulator.h +++ b/tests/lbpm_color_simulator.h @@ -38,40 +38,54 @@ private: // Helper class to compute the blob ids +static const std::string id_map_filename = "lbpm_id_map.txt"; typedef std::shared_ptr > BlobIDstruct; +typedef std::shared_ptr > BlobIDList; class BlobIdentificationWorkItem: public ThreadPool::WorkItem { public: - BlobIdentificationWorkItem( int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, + BlobIdentificationWorkItem( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, std::shared_ptr phase_, const DoubleArray& dist_, - BlobIDstruct last_id_, BlobIDstruct new_id_ ): - Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), phase(phase_), - dist(dist_), last_id(last_id_), new_id(new_id_) { } + BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ): + timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), + phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_) { } virtual void run() { ThreadPool::WorkItem::d_state = 1; // Change state to in progress // Compute the global blob id and compare to the previous version PROFILE_START("Identify blobs and maps",1); double vF = 0.0; double vS = 0.0; - IntArray& ids = new_id->second; - new_id->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids); - if ( last_id==NULL ) { + IntArray& ids = new_index->second; + new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids); + static int max_id = -1; + new_id->first = new_index->first; + new_id->second = new_index->second; + if ( last_id!=NULL ) { // Compute the timestep-timestep map - const IntArray& old_ids = new_id->second; + const IntArray& old_ids = last_id->second; ID_map_struct map = computeIDMap(old_ids,ids); // Renumber the current timestep's ids - + getNewIDs(map,max_id,*new_list); + renumberIDs(*new_list,new_id->second); + writeIDMap(map,timestep,id_map_filename); + } else { + max_id = -1; + ID_map_struct map(new_id->first); + getNewIDs(map,max_id,*new_list); + writeIDMap(map,timestep,id_map_filename); } PROFILE_STOP("Identify blobs and maps",1); ThreadPool::WorkItem::d_state = 2; // Change state to finished } private: BlobIdentificationWorkItem(); + int timestep; int Nx, Ny, Nz; const RankInfoStruct& rank_info; std::shared_ptr phase; const DoubleArray& dist; - BlobIDstruct last_id, new_id; + BlobIDstruct last_id, new_index, new_id; + BlobIDList new_list; }; @@ -80,12 +94,15 @@ private: class AnalysisWorkItem: public ThreadPool::WorkItem { public: - AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, BlobIDstruct ids, double beta_ ): - type(type_), timestep(timestep_), Averages(Averages_), blob_ids(ids), beta(beta_) { } + AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, + BlobIDstruct ids, BlobIDList id_list_, double beta_ ): + type(type_), timestep(timestep_), Averages(Averages_), + blob_ids(ids), id_list(id_list_), beta(beta_) { } virtual void run() { ThreadPool::WorkItem::d_state = 1; // Change state to in progress Averages.NumberComponents_NWP = blob_ids->first; Averages.Label_NWP = blob_ids->second; + Averages.Label_NWP_map = *id_list; Averages.NumberComponents_WP = 1; Averages.Label_WP.fill(0.0); if ( (type&CopyPhaseIndicator) != 0 ) { @@ -115,13 +132,15 @@ private: int timestep; TwoPhase& Averages; BlobIDstruct blob_ids; + BlobIDList id_list; double beta; }; // Function to start the analysis void run_analysis( int timestep, int restart_interval, - const RankInfoStruct& rank_info, TwoPhase& Averages, BlobIDstruct& last_ids, + const RankInfoStruct& rank_info, TwoPhase& Averages, + BlobIDstruct& last_ids, BlobIDstruct& last_index, BlobIDList& last_id_map, int Nx, int Ny, int Nz, bool pBC, double beta, double err, const double *Phi, double *Pressure, const double *Velocity, const char *ID, const double *f_even, const double *f_odd, const double *Den, @@ -191,17 +210,22 @@ void run_analysis( int timestep, int restart_interval, // Spawn threads to do blob identification work if ( (type&IdentifyBlobs)!=0 ) { + BlobIDstruct new_index(new std::pair(0,IntArray())); BlobIDstruct new_ids(new std::pair(0,IntArray())); - ThreadPool::WorkItem *work = new BlobIdentificationWorkItem( - Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_ids); + BlobIDList new_list(new std::vector()); + ThreadPool::WorkItem *work = new BlobIdentificationWorkItem(timestep, + Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_index,new_ids,new_list); work->add_dependency(wait.blobID); + last_index = new_index; last_ids = new_ids; + last_id_map = new_list; wait.blobID = tpool.add_work(work); } // Spawn threads to do the analysis work if ( (type&CalcDist) != 0 ) { - ThreadPool::WorkItem *work = new AnalysisWorkItem(type,timestep,Averages,last_ids,beta); + ThreadPool::WorkItem *work = new AnalysisWorkItem( + type,timestep,Averages,last_index,last_id_map,beta); work->add_dependency(wait.blobID); work->add_dependency(wait.analysis); wait.analysis = tpool.add_work(work); diff --git a/visit/avtLBMFileFormat.C b/visit/avtLBMFileFormat.C index 758991bb..2cf03955 100644 --- a/visit/avtLBMFileFormat.C +++ b/visit/avtLBMFileFormat.C @@ -99,9 +99,23 @@ #endif +// Null stream buffer +class NullBufferClass : public std::streambuf { +public: + virtual std::streamsize xsputn(const char *s, std::streamsize n) { return n; } + virtual int overflow(int c) { return 1; } + virtual void setOutputStream( std::ostream* ) {} +}; + + // Output streams -static IO::ParallelStreamBuffer DebugStreamBuffer1; -static IO::ParallelStreamBuffer DebugStreamBuffer2; +#if 0 + static IO::ParallelStreamBuffer DebugStreamBuffer1; + static IO::ParallelStreamBuffer DebugStreamBuffer2; +#else + static NullBufferClass DebugStreamBuffer1; + static NullBufferClass DebugStreamBuffer2; +#endif std::ostream DebugStream1(&DebugStreamBuffer1); std::ostream DebugStream2(&DebugStreamBuffer2); diff --git a/visit/vtkHelpers.h b/visit/vtkHelpers.h index 04d347c7..ba1e8b85 100644 --- a/visit/vtkHelpers.h +++ b/visit/vtkHelpers.h @@ -196,7 +196,7 @@ inline vtkDataSet* meshToVTK( std::shared_ptr mesh ) mesh2 = DomainToVTK( std::dynamic_pointer_cast(mesh) ); DebugStream1 << " Volume mesh created" << std::endl; } else { - DebugStream1 << " Error, unknown mesh type" << std::endl; + //DebugStream1 << " Error, unknown mesh type" << std::endl; return NULL; } return mesh2;