diff --git a/.clang-format b/.clang-format new file mode 100644 index 00000000..b930c24c --- /dev/null +++ b/.clang-format @@ -0,0 +1,108 @@ +# To run clang tools: +# cd to root directory +# To update format only: +# find . -name "*.cpp" -or -name "*.cc" -or -name "*.h" -or -name "*.hpp" -or -name "*.I" | xargs -I{} clang-format -i {} +# git status -s . | sed s/^...// | grep -E "(\.cpp|\.h|\.cc|\.hpp|\.I)" | xargs -I{} clang-format -i {} + +# To run modernize +# export CLANG_PATH=/packages/llvm/build/llvm-60 +# export PATH=${CLANG_PATH}/bin:${CLANG_PATH}/share/clang:$PATH +# find src -name "*.cpp" -or -name "*.cc" | xargs -I{} clang-tidy -checks=modernize* -p=/projects/AtomicModel/build/debug -fix {} +# find src -name "*.cpp" -or -name "*.cc" -or -name "*.h" -or -name "*.hpp" -or -name "*.I" | xargs -I{} clang-format -i {} + + + + + +--- +Language: Cpp +# BasedOnStyle: LLVM +AccessModifierOffset: -4 +AlignAfterOpenBracket: DontAlign +AlignConsecutiveAssignments: true +AlignConsecutiveDeclarations: false +AlignEscapedNewlinesLeft: true +AlignOperands: true +AlignTrailingComments: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortBlocksOnASingleLine: true +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: false +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: true +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterClass: true + AfterControlStatement: false + AfterEnum: false + AfterFunction: true + AfterNamespace: false + AfterObjCDeclaration: true + AfterStruct: false + AfterUnion: false + BeforeCatch: false + BeforeElse: false + IndentBraces: false +BreakBeforeBinaryOperators: None +#BreakBeforeBraces: Stroustrup +BreakBeforeBraces: Custom +BreakBeforeTernaryOperators: false +BreakConstructorInitializersBeforeComma: false +ColumnLimit: 100 +CommentPragmas: '^ IWYU pragma:' +ConstructorInitializerAllOnOneLineOrOnePerLine: true +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: false +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: false +ForEachMacros: [ foreach, Q_FOREACH, BOOST_FOREACH ] +IncludeCategories: + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 2 + - Regex: '^(<|"(gtest|isl|json)/)' + Priority: 3 + - Regex: '.*' + Priority: 1 +IndentCaseLabels: false +IndentWidth: 4 +IndentWrappedFunctionNames: false +KeepEmptyLinesAtTheStartOfBlocks: true +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 2 +NamespaceIndentation: None +ObjCBlockIndentWidth: 4 +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyExcessCharacter: 1000 +PenaltyReturnTypeOnItsOwnLine: 60 +PointerAlignment: Right +ReflowComments: true +SortIncludes: true +SortUsingDeclarations: true +SpaceAfterCStyleCast: true +SpaceAfterTemplateKeyword: false +SpaceBeforeAssignmentOperators: true +SpaceBeforeParens: ControlStatements +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInParentheses: true +SpacesInSquareBrackets: false +Standard: Cpp11 +TabWidth: 4 +UseTab: Never +... + diff --git a/IO/silo.hpp b/IO/silo.hpp index 35852004..1e17aa5c 100644 --- a/IO/silo.hpp +++ b/IO/silo.hpp @@ -2,8 +2,8 @@ #define SILO_INTERFACE_HPP #include "IO/silo.h" -#include "common/Utilities.h" #include "common/MPI.h" +#include "common/Utilities.h" #include "ProfilerApp.h" @@ -13,52 +13,77 @@ #include - namespace silo { /**************************************************** -* Helper functions * -****************************************************/ -template static constexpr int getType(); -template<> constexpr int getType() { return DB_DOUBLE; } -template<> constexpr int getType() { return DB_FLOAT; } -template<> constexpr int getType() { return DB_INT; } + * Helper functions * + ****************************************************/ template -inline void copyData( Array& data, int type, const void *src ) +static constexpr int getType(); +template<> +constexpr int getType() +{ + return DB_DOUBLE; +} +template<> +constexpr int getType() +{ + return DB_FLOAT; +} +template<> +constexpr int getType() +{ + return DB_INT; +} +template +inline void copyData( Array &data, int type, const void *src ) { if ( type == getType() ) - memcpy( data.data(), src, data.length()*sizeof(TYPE) ); + memcpy( data.data(), src, data.length() * sizeof( TYPE ) ); else if ( type == DB_DOUBLE ) - data.copy( static_cast(src) ); + data.copy( static_cast( src ) ); else if ( type == DB_FLOAT ) - data.copy( static_cast(src) ); + data.copy( static_cast( src ) ); else if ( type == DB_INT ) - data.copy( static_cast(src) ); + data.copy( static_cast( src ) ); else - ERROR("Unknown type"); + ERROR( "Unknown type" ); } /**************************************************** -* Write/read an arbitrary vector * -****************************************************/ -template constexpr int getSiloType(); -template<> constexpr int getSiloType() { return DB_INT; } -template<> constexpr int getSiloType() { return DB_FLOAT; } -template<> constexpr int getSiloType() { return DB_DOUBLE; } + * Write/read an arbitrary vector * + ****************************************************/ template -void write( DBfile* fid, const std::string& varname, const std::vector& data ) +constexpr int getSiloType(); +template<> +constexpr int getSiloType() +{ + return DB_INT; +} +template<> +constexpr int getSiloType() +{ + return DB_FLOAT; +} +template<> +constexpr int getSiloType() +{ + return DB_DOUBLE; +} +template +void write( DBfile *fid, const std::string &varname, const std::vector &data ) { int dims = data.size(); - int err = DBWrite( fid, varname.c_str(), (void*) data.data(), &dims, 1, getSiloType() ); + int err = DBWrite( fid, varname.c_str(), (void *) data.data(), &dims, 1, getSiloType() ); ASSERT( err == 0 ); } template -std::vector read( DBfile* fid, const std::string& varname ) +std::vector read( DBfile *fid, const std::string &varname ) { int N = DBGetVarLength( fid, varname.c_str() ); - std::vector data(N); + std::vector data( N ); int err = DBReadVar( fid, varname.c_str(), data.data() ); ASSERT( err == 0 ); return data; @@ -66,31 +91,31 @@ std::vector read( DBfile* fid, const std::string& varname ) /**************************************************** -* Helper function to get variable suffixes * -****************************************************/ + * Helper function to get variable suffixes * + ****************************************************/ inline std::vector getVarSuffix( int ndim, int nvars ) { - std::vector suffix(nvars); + std::vector suffix( nvars ); if ( nvars == 1 ) { suffix[0] = ""; } else if ( nvars == ndim ) { - if ( ndim==2 ) { + if ( ndim == 2 ) { suffix[0] = "_x"; suffix[1] = "_y"; - } else if ( ndim==3 ) { + } else if ( ndim == 3 ) { suffix[0] = "_x"; suffix[1] = "_y"; suffix[2] = "_z"; } else { - ERROR("Not finished"); + ERROR( "Not finished" ); } - } else if ( nvars == ndim*ndim ) { - if ( ndim==2 ) { + } else if ( nvars == ndim * ndim ) { + if ( ndim == 2 ) { suffix[0] = "_xx"; suffix[1] = "_xy"; suffix[2] = "_yx"; suffix[3] = "_yy"; - } else if ( ndim==3 ) { + } else if ( ndim == 3 ) { suffix[0] = "_xx"; suffix[1] = "_xy"; suffix[2] = "_xz"; @@ -101,122 +126,127 @@ inline std::vector getVarSuffix( int ndim, int nvars ) suffix[7] = "_zy"; suffix[8] = "_zz"; } else { - ERROR("Not finished"); + ERROR( "Not finished" ); } } else { - for (int i=0; i -void writeUniformMesh( DBfile* fid, const std::string& meshname, - const std::array& range, const std::array& N ) +void writeUniformMesh( DBfile *fid, const std::string &meshname, + const std::array &range, const std::array &N ) { - PROFILE_START("writeUniformMesh",2); + PROFILE_START( "writeUniformMesh", 2 ); int dims[NDIM]; - for (size_t d=0; d= 1 ) { x = new float[dims[0]]; - for (int i=0; i= 2 ) { y = new float[dims[1]]; - for (int i=0; i= 3 ) { z = new float[dims[2]]; - for (int i=0; i -void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const std::array& N, - const std::string& varname, const Array& data, VariableType type ) + * Write a vector/tensor quad variable * + ****************************************************/ +template +void writeUniformMeshVariable( DBfile *fid, const std::string &meshname, + const std::array &N, const std::string &varname, const Array &data, + VariableType type ) { - PROFILE_START("writeUniformMeshVariable",2); - int nvars=1, dims[NDIM]={1}; + PROFILE_START( "writeUniformMeshVariable", 2 ); + int nvars = 1, dims[NDIM] = { 1 }; const TYPE *vars[NDIM] = { nullptr }; - int vartype = 0; + int vartype = 0; if ( type == VariableType::NodeVariable ) { - ASSERT( data.ndim()==NDIM || data.ndim()==NDIM+1 ); - for (int d=0; d var_names(nvars); - for (int i=0; i var_names( nvars ); + for ( int i = 0; i < nvars; i++ ) var_names[i] = varname + suffix[i]; - std::vector varnames(nvars,nullptr); - for (int i=0; i(var_names[i].c_str()); - int err = DBPutQuadvar( fid, varname.c_str(), meshname.c_str(), nvars, - varnames.data(), vars, dims, NDIM, nullptr, 0, getType(), vartype, nullptr ); + std::vector varnames( nvars, nullptr ); + for ( int i = 0; i < nvars; i++ ) + varnames[i] = const_cast( var_names[i].c_str() ); + int err = DBPutQuadvar( fid, varname.c_str(), meshname.c_str(), nvars, varnames.data(), vars, + dims, NDIM, nullptr, 0, getType(), vartype, nullptr ); ASSERT( err == 0 ); - PROFILE_STOP("writeUniformMeshVariable",2); + PROFILE_STOP( "writeUniformMeshVariable", 2 ); } -template -Array readUniformMeshVariable( DBfile* fid, const std::string& varname ) +template +Array readUniformMeshVariable( DBfile *fid, const std::string &varname ) { auto var = DBGetQuadvar( fid, varname.c_str() ); ASSERT( var != nullptr ); Array data( var->nels, var->nvals ); int type = var->datatype; - for (int i=0; invals; i++) { + for ( int i = 0; i < var->nvals; i++ ) { Array data2( var->nels ); copyData( data2, type, var->vals[i] ); - memcpy( &data(0,i), data2.data(), var->nels*sizeof(TYPE) ); + memcpy( &data( 0, i ), data2.data(), var->nels * sizeof( TYPE ) ); } DBFreeQuadvar( var ); - std::vector dims( var->ndims+1, var->nvals ); - for (int d=0; dndims; d++) + std::vector dims( var->ndims + 1, var->nvals ); + for ( int d = 0; d < var->ndims; d++ ) dims[d] = var->dims[d]; data.reshape( dims ); return data; @@ -224,54 +254,55 @@ Array readUniformMeshVariable( DBfile* fid, const std::string& varname ) /**************************************************** -* Read/write a point mesh/variable to silo * -****************************************************/ + * Read/write a point mesh/variable to silo * + ****************************************************/ template -void writePointMesh( DBfile* fid, const std::string& meshname, - int ndim, int N, const TYPE *coords[] ) +void writePointMesh( + DBfile *fid, const std::string &meshname, int ndim, int N, const TYPE *coords[] ) { int err = DBPutPointmesh( fid, meshname.c_str(), ndim, coords, N, getType(), nullptr ); ASSERT( err == 0 ); } -template -Array readPointMesh( DBfile* fid, const std::string& meshname ) +template +Array readPointMesh( DBfile *fid, const std::string &meshname ) { auto mesh = DBGetPointmesh( fid, meshname.c_str() ); - int N = mesh->nels; - int ndim = mesh->ndims; - Array coords(N,ndim); + int N = mesh->nels; + int ndim = mesh->ndims; + Array coords( N, ndim ); int type = mesh->datatype; - for (int d=0; d data2( N ); copyData( data2, type, mesh->coords[d] ); - memcpy( &coords(0,d), data2.data(), N*sizeof(TYPE) ); + memcpy( &coords( 0, d ), data2.data(), N * sizeof( TYPE ) ); } DBFreePointmesh( mesh ); return coords; } template -void writePointMeshVariable( DBfile* fid, const std::string& meshname, - const std::string& varname, const Array& data ) +void writePointMeshVariable( + DBfile *fid, const std::string &meshname, const std::string &varname, const Array &data ) { - int N = data.size(0); - int nvars = data.size(1); - std::vector vars(nvars); - for (int i=0; i(), nullptr ); + int N = data.size( 0 ); + int nvars = data.size( 1 ); + std::vector vars( nvars ); + for ( int i = 0; i < nvars; i++ ) + vars[i] = &data( 0, i ); + int err = DBPutPointvar( + fid, varname.c_str(), meshname.c_str(), nvars, vars.data(), N, getType(), nullptr ); ASSERT( err == 0 ); } -template -Array readPointMeshVariable( DBfile* fid, const std::string& varname ) +template +Array readPointMeshVariable( DBfile *fid, const std::string &varname ) { auto var = DBGetPointvar( fid, varname.c_str() ); ASSERT( var != nullptr ); Array data( var->nels, var->nvals ); int type = var->datatype; - for (int i=0; invals; i++) { + for ( int i = 0; i < var->nvals; i++ ) { Array data2( var->nels ); copyData( data2, type, var->vals[i] ); - memcpy( &data(0,i), data2.data(), var->nels*sizeof(TYPE) ); + memcpy( &data( 0, i ), data2.data(), var->nels * sizeof( TYPE ) ); } DBFreeMeshvar( var ); return data; @@ -279,110 +310,110 @@ Array readPointMeshVariable( DBfile* fid, const std::string& varname ) /**************************************************** -* Read/write a triangle mesh * -****************************************************/ + * Read/write a triangle mesh * + ****************************************************/ template -void writeTriMesh( DBfile* fid, const std::string& meshName, - int ndim, int ndim_tri, int N, const TYPE *coords[], int N_tri, const int *tri[] ) +void writeTriMesh( DBfile *fid, const std::string &meshName, int ndim, int ndim_tri, int N, + const TYPE *coords[], int N_tri, const int *tri[] ) { auto zoneName = meshName + "_zones"; - std::vector nodelist( (ndim_tri+1)*N_tri ); - for (int i=0, j=0; i nodelist( ( ndim_tri + 1 ) * N_tri ); + for ( int i = 0, j = 0; i < N_tri; i++ ) { + for ( int d = 0; d < ndim_tri + 1; d++, j++ ) nodelist[j] = tri[d][i]; } int shapetype = 0; - if ( ndim_tri==1 ) + if ( ndim_tri == 1 ) shapetype = DB_ZONETYPE_BEAM; - else if ( ndim_tri==2 ) + else if ( ndim_tri == 2 ) shapetype = DB_ZONETYPE_TRIANGLE; - else if ( ndim_tri==3 ) + else if ( ndim_tri == 3 ) shapetype = DB_ZONETYPE_PYRAMID; else - ERROR("Unknown shapetype"); - int shapesize = ndim_tri+1; - int shapecnt = N_tri; - DBPutZonelist2( fid, zoneName.c_str(), N_tri, ndim_tri, nodelist.data(), - nodelist.size(), 0, 0, 0, &shapetype, &shapesize, &shapecnt, 1, nullptr ); - DBPutUcdmesh( fid, meshName.c_str(), ndim, nullptr, coords, N, - nodelist.size(), zoneName.c_str(), nullptr, getType(), nullptr ); + ERROR( "Unknown shapetype" ); + int shapesize = ndim_tri + 1; + int shapecnt = N_tri; + DBPutZonelist2( fid, zoneName.c_str(), N_tri, ndim_tri, nodelist.data(), nodelist.size(), 0, 0, + 0, &shapetype, &shapesize, &shapecnt, 1, nullptr ); + DBPutUcdmesh( fid, meshName.c_str(), ndim, nullptr, coords, N, nodelist.size(), + zoneName.c_str(), nullptr, getType(), nullptr ); } template -void readTriMesh( DBfile* fid, const std::string& meshname, Array& coords, Array& tri ) +void readTriMesh( DBfile *fid, const std::string &meshname, Array &coords, Array &tri ) { - auto mesh = DBGetUcdmesh( fid, meshname.c_str() ); - int ndim = mesh->ndims; + auto mesh = DBGetUcdmesh( fid, meshname.c_str() ); + int ndim = mesh->ndims; int N_nodes = mesh->nnodes; - coords.resize(N_nodes,ndim); + coords.resize( N_nodes, ndim ); int mesh_type = mesh->datatype; - for (int d=0; d data2( N_nodes ); copyData( data2, mesh_type, mesh->coords[d] ); - memcpy( &coords(0,d), data2.data(), N_nodes*sizeof(TYPE) ); + memcpy( &coords( 0, d ), data2.data(), N_nodes * sizeof( TYPE ) ); } - auto zones = mesh->zones; + auto zones = mesh->zones; int N_zones = zones->nzones; - ASSERT( zones->nshapes==1 ); + ASSERT( zones->nshapes == 1 ); int shapesize = zones->shapesize[0]; - tri.resize(N_zones,shapesize); - for (int i=0; inodelist[i*shapesize+j]; + tri.resize( N_zones, shapesize ); + for ( int i = 0; i < N_zones; i++ ) { + for ( int j = 0; j < shapesize; j++ ) + tri( i, j ) = zones->nodelist[i * shapesize + j]; } DBFreeUcdmesh( mesh ); } template -void writeTriMeshVariable( DBfile* fid, int ndim, const std::string& meshname, - const std::string& varname, const Array& data, VariableType type ) +void writeTriMeshVariable( DBfile *fid, int ndim, const std::string &meshname, + const std::string &varname, const Array &data, VariableType type ) { - int nvars = 0; - int vartype = 0; + int nvars = 0; + int vartype = 0; const TYPE *vars[10] = { nullptr }; if ( type == VariableType::NodeVariable ) { vartype = DB_NODECENT; - nvars = data.size(1); - for (int i=0; i var_names(nvars); - for (int i=0; i var_names( nvars ); + for ( int i = 0; i < nvars; i++ ) var_names[i] = varname + suffix[i]; - std::vector varnames(nvars,nullptr); - for (int i=0; i(var_names[i].c_str()); - DBPutUcdvar( fid, varname.c_str(), meshname.c_str(), nvars, - varnames.data(), vars, data.size(0), nullptr, 0, getType(), vartype, nullptr ); + std::vector varnames( nvars, nullptr ); + for ( int i = 0; i < nvars; i++ ) + varnames[i] = const_cast( var_names[i].c_str() ); + DBPutUcdvar( fid, varname.c_str(), meshname.c_str(), nvars, varnames.data(), vars, + data.size( 0 ), nullptr, 0, getType(), vartype, nullptr ); } template -Array readTriMeshVariable( DBfile* fid, const std::string& varname ) +Array readTriMeshVariable( DBfile *fid, const std::string &varname ) { auto var = DBGetUcdvar( fid, varname.c_str() ); ASSERT( var != nullptr ); Array data( var->nels, var->nvals ); int type = var->datatype; - for (int i=0; invals; i++) { + for ( int i = 0; i < var->nvals; i++ ) { Array data2( var->nels ); copyData( data2, type, var->vals[i] ); - memcpy( &data(0,i), data2.data(), var->nels*sizeof(TYPE) ); + memcpy( &data( 0, i ), data2.data(), var->nels * sizeof( TYPE ) ); } DBFreeUcdvar( var ); return data; } -}; // silo namespace +}; // namespace silo #endif diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 79ef5c93..d0657391 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -156,7 +156,7 @@ void SubPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB, dou } void SubPhase::Basic(){ - int i,j,k,n,imin,jmin,kmin,kmax, nq; + int i,j,k,n,imin,jmin,kmin,kmax; // If external boundary conditions are set, do not average over the inlet kmin=1; kmax=Nz-1; @@ -172,8 +172,20 @@ void SubPhase::Basic(){ double count_w = 0.0; double count_n = 0.0; - total_wetting_interaction = count_wetting_interaction = 0.0; - total_wetting_interaction_global = count_wetting_interaction_global=0.0; + /* compute the laplacian */ + Dm->CommunicateMeshHalo(Phi); + for (int k=1; kCommunicateMeshHalo(DelPhi); for (k=0; kid[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 1\n",i,j,k); - local_wetting_interaction += (Phi(nq)-wetval); - local_wetting_weight += 1.0; - } - nq = (k)*Nx*Ny+(j)*Nx+(i-1); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 2\n",i,j,k); - local_wetting_interaction += (Phi(nq)-wetval); - local_wetting_weight += 1.0; - } - nq = (k)*Nx*Ny+(j+1)*Nx+(i); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 3\n",i,j,k); - local_wetting_interaction += (Phi(nq)-wetval); - local_wetting_weight += 1.0; - } - nq = (k)*Nx*Ny+(j-1)*Nx+(i); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 4\n",i,j,k); - local_wetting_interaction += (Phi(nq)-wetval); - local_wetting_weight += 1.0; - } - nq = (k+1)*Nx*Ny+(j)*Nx+(i); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 5\n",i,j,k); - local_wetting_interaction += (Phi(nq)-wetval); - local_wetting_weight += 1.0; - } - nq = (k-1)*Nx*Ny+(j)*Nx+(i); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 6\n",i,j,k); - local_wetting_interaction += (Phi(nq)-wetval); - local_wetting_weight += 1.0; - } - // x, y interactions - nq = (k)*Nx*Ny+(j+1)*Nx+(i+1); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 7\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - nq = (k)*Nx*Ny+(j-1)*Nx+(i-1); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 8\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - nq = (k)*Nx*Ny+(j-1)*Nx+(i+1); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 9\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - nq = (k)*Nx*Ny+(j+1)*Nx+(i-1); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 10\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - // xz interactions - nq = (k+1)*Nx*Ny+(j)*Nx+(i+1); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 11\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - nq = (k-1)*Nx*Ny+(j)*Nx+(i-1); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 12\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - nq = (k+1)*Nx*Ny+(j)*Nx+(i-1); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 13\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - nq = (k-1)*Nx*Ny+(j)*Nx+(i+1); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 14\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - // yz interactions - nq = (k+1)*Nx*Ny+(j+1)*Nx+(i); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 15\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - nq = (k-1)*Nx*Ny+(j-1)*Nx+(i); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 16\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - nq = (k+1)*Nx*Ny+(j-1)*Nx+(i); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 17\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - nq = (k-1)*Nx*Ny+(j+1)*Nx+(i); - if ( Dm->id[nq] > 0 ) { - //if (Phi(nq)!=Phi(nq)) printf("%i %i %i : 18\n",i,j,k); - local_wetting_interaction += 0.5*(Phi(nq)-wetval); - local_wetting_weight += 0.5; - } - /* interaction due to this solid site*/ - if (local_wetting_interaction == local_wetting_interaction){ - total_wetting_interaction += 0.5*local_wetting_interaction; - if (local_wetting_weight > 0.0) - count_wetting_interaction += local_wetting_weight; - } - else{ - //printf("Check interaction at %i %i %i \n",i,j,k); - } + } + } + } + + total_wetting_interaction = count_wetting_interaction = 0.0; + total_wetting_interaction_global = count_wetting_interaction_global=0.0; + for (k=kmin; kid[n] > 0 && SDs(i,j,k) < 2.0 ){ + count_wetting_interaction += 1.0; + total_wetting_interaction += DelPhi(i,j,k); } } } @@ -367,11 +267,12 @@ void SubPhase::Basic(){ //printf("wetting interaction = %f, count = %f\n",total_wetting_interaction,count_wetting_interaction); total_wetting_interaction_global=Dm->Comm.sumReduce( total_wetting_interaction); count_wetting_interaction_global=Dm->Comm.sumReduce( count_wetting_interaction); - /* normalize wetting interactions */ + /* normalize wetting interactions <-- Don't do this if normalizing laplacian (use solid surface area) if (count_wetting_interaction > 0.0) total_wetting_interaction /= count_wetting_interaction; if (count_wetting_interaction_global > 0.0) total_wetting_interaction_global /= count_wetting_interaction_global; + */ gwb.V=Dm->Comm.sumReduce( wb.V); gnb.V=Dm->Comm.sumReduce( nb.V); diff --git a/analysis/SubPhase.h b/analysis/SubPhase.h index d2ce6b44..a6d35edd 100644 --- a/analysis/SubPhase.h +++ b/analysis/SubPhase.h @@ -68,11 +68,11 @@ public: * b - bulk (total) */ // local entities - phase wc,wd,wb,nc,nd,nb; + phase wc,wd,wb,nc,nd,nb,solid; interface iwn,iwnc; // global entities - phase gwc,gwd,gwb,gnc,gnd,gnb; + phase gwc,gwd,gwb,gnc,gnd,gnb,gsolid; interface giwn,giwnc; /* fluid-solid wetting interaction */ double total_wetting_interaction, count_wetting_interaction; diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index ad231f3f..37f58d0c 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -542,7 +542,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptrrank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); Dm->Comm.barrier(); diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index e2fca48d..f43a26ff 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -13,140 +13,170 @@ #include "ProfilerApp.h" -AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs) +AnalysisType &operator|=( AnalysisType &lhs, AnalysisType rhs ) { - lhs = static_cast( - static_cast::type>(lhs) | - static_cast::type>(rhs) - ); + lhs = static_cast( static_cast::type>( lhs ) | + static_cast::type>( rhs ) ); return lhs; } bool matches( AnalysisType x, AnalysisType y ) { - return ( static_cast::type>(x) & - static_cast::type>(y) ) != 0; + return ( static_cast::type>( x ) & + static_cast::type>( y ) ) != 0; } +// Create a shared_ptr to an array of values template -void DeleteArray( const TYPE *p ) +static inline std::shared_ptr make_shared_array( size_t N ) { - delete [] p; + return std::shared_ptr( new TYPE[N], []( const TYPE *p ) { delete[] p; } ); } // Helper class to write the restart file from a seperate thread -class WriteRestartWorkItem: public ThreadPool::WorkItemRet +class WriteRestartWorkItem : public ThreadPool::WorkItemRet { public: - WriteRestartWorkItem( const char* filename_, std::shared_ptr cDen_, std::shared_ptr cfq_, int N_ ): - filename(filename_), cfq(cfq_), cDen(cDen_), N(N_) {} - virtual void run() { - PROFILE_START("Save Checkpoint",1); + WriteRestartWorkItem( const std::string &filename_, std::shared_ptr cDen_, + std::shared_ptr cfq_, int N_ ) + : filename( filename_ ), cfq( cfq_ ), cDen( cDen_ ), N( N_ ) + { + } + virtual void run() + { + PROFILE_START( "Save Checkpoint", 1 ); double value; - ofstream File(filename,ios::binary); - for (int n=0; n cfq,cDen; - // const DoubleArray& phase; - //const DoubleArray& dist; + const std::string filename; + std::shared_ptr cfq, cDen; const int N; }; // Helper class to compute the blob ids +typedef std::shared_ptr> BlobIDstruct; +typedef std::shared_ptr> BlobIDList; static const std::string id_map_filename = "lbpm_id_map.txt"; -class BlobIdentificationWorkItem1: public ThreadPool::WorkItemRet +class BlobIdentificationWorkItem1 : public ThreadPool::WorkItemRet { public: - BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, - std::shared_ptr phase_, const DoubleArray& dist_, - BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper&& comm_ ): - 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_), comm(std::move(comm_)) -{ -} - ~BlobIdentificationWorkItem1() { } - virtual void run() { - // Compute the global blob id and compare to the previous version - PROFILE_START("Identify blobs",1); - double vF = 0.0; - double vS = -1.0; // one voxel buffer region around solid - IntArray& ids = new_index->second; - new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,comm.comm); - PROFILE_STOP("Identify blobs",1); + BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, + const RankInfoStruct &rank_info_, std::shared_ptr phase_, + const DoubleArray &dist_, BlobIDstruct last_id_, BlobIDstruct new_index_, + BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper &&comm_ ) + : 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_ ), + comm( std::move( comm_ ) ) + { } + ~BlobIdentificationWorkItem1() {} + virtual void run() + { + // Compute the global blob id and compare to the previous version + PROFILE_START( "Identify blobs", 1 ); + double vF = 0.0; + double vS = -1.0; // one voxel buffer region around solid + IntArray &ids = new_index->second; + new_index->first = ComputeGlobalBlobIDs( + Nx - 2, Ny - 2, Nz - 2, rank_info, *phase, dist, vF, vS, ids, comm.comm ); + PROFILE_STOP( "Identify blobs", 1 ); + } + private: BlobIdentificationWorkItem1(); int timestep; int Nx, Ny, Nz; - const RankInfoStruct& rank_info; + const RankInfoStruct rank_info; std::shared_ptr phase; - const DoubleArray& dist; + const DoubleArray &dist; BlobIDstruct last_id, new_index, new_id; BlobIDList new_list; runAnalysis::commWrapper comm; }; -class BlobIdentificationWorkItem2: public ThreadPool::WorkItemRet +class BlobIdentificationWorkItem2 : public ThreadPool::WorkItemRet { public: - BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, - std::shared_ptr phase_, const DoubleArray& dist_, - BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ , runAnalysis::commWrapper&& comm_ ): - 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_), comm(std::move(comm_)) -{ -} - ~BlobIdentificationWorkItem2() { } - virtual void run() { + BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, + const RankInfoStruct &rank_info_, std::shared_ptr phase_, + const DoubleArray &dist_, BlobIDstruct last_id_, BlobIDstruct new_index_, + BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper &&comm_ ) + : 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_ ), + comm( std::move( comm_ ) ) + { + } + ~BlobIdentificationWorkItem2() {} + virtual void run() + { // Compute the global blob id and compare to the previous version - PROFILE_START("Identify blobs maps",1); - const IntArray& ids = new_index->second; - static int max_id = -1; - new_id->first = new_index->first; - new_id->second = new_index->second; - if ( last_id.get()!=NULL ) { + PROFILE_START( "Identify blobs maps", 1 ); + const IntArray &ids = new_index->second; + static int max_id = -1; + new_id->first = new_index->first; + new_id->second = new_index->second; + if ( last_id.get() != NULL ) { // Compute the timestep-timestep map - const IntArray& old_ids = last_id->second; - ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,comm.comm); + const IntArray &old_ids = last_id->second; + ID_map_struct map = computeIDMap( Nx, Ny, Nz, old_ids, ids, comm.comm ); // Renumber the current timestep's ids - getNewIDs(map,max_id,*new_list); - renumberIDs(*new_list,new_id->second); - writeIDMap(map,timestep,id_map_filename); + 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); + ID_map_struct map( new_id->first ); + getNewIDs( map, max_id, *new_list ); + writeIDMap( map, timestep, id_map_filename ); } - PROFILE_STOP("Identify blobs maps",1); + PROFILE_STOP( "Identify blobs maps", 1 ); } + private: BlobIdentificationWorkItem2(); int timestep; int Nx, Ny, Nz; - const RankInfoStruct& rank_info; + const RankInfoStruct rank_info; std::shared_ptr phase; - const DoubleArray& dist; + const DoubleArray &dist; BlobIDstruct last_id, new_index, new_id; BlobIDList new_list; runAnalysis::commWrapper comm; @@ -154,324 +184,375 @@ private: // Helper class to write the vis file from a thread -class WriteVisWorkItem: public ThreadPool::WorkItemRet +class WriteVisWorkItem : public ThreadPool::WorkItemRet { public: - WriteVisWorkItem( int timestep_, std::vector& visData_, - TwoPhase& Avgerages_, fillHalo& fillData_, runAnalysis::commWrapper&& comm_ ): - timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_), comm(std::move(comm_)) - { - } - ~WriteVisWorkItem() { } - virtual void run() { - PROFILE_START("Save Vis",1); - - ASSERT(visData[0].vars[0]->name=="phase"); - Array& PhaseData = visData[0].vars[0]->data; - fillData.copy(Averages.SDn,PhaseData); + WriteVisWorkItem( int timestep_, std::vector &visData_, + TwoPhase &Avgerages_, std::array n_, RankInfoStruct rank_info_, + runAnalysis::commWrapper &&comm_ ) + : timestep( timestep_ ), + visData( visData_ ), + Averages( Avgerages_ ), + n( std::move( n_ ) ), + rank_info( std::move( rank_info_ ) ), + comm( std::move( comm_ ) ) + { + } + ~WriteVisWorkItem() {} + virtual void run() + { + PROFILE_START( "Save Vis", 1 ); - ASSERT(visData[0].vars[5]->name=="SignDist"); - Array& SignData = visData[0].vars[5]->data; - fillData.copy(Averages.SDs,SignData); + fillHalo fillData( comm.comm, rank_info, n, { 1, 1, 1 }, 0, 1 ); - ASSERT(visData[0].vars[1]->name=="Pressure"); - Array& PressData = visData[0].vars[1]->data; - fillData.copy(Averages.Press,PressData); + ASSERT( visData[0].vars[0]->name == "phase" ); + Array &PhaseData = visData[0].vars[0]->data; + fillData.copy( Averages.SDn, PhaseData ); - ASSERT(visData[0].vars[2]->name=="Velocity_x"); - ASSERT(visData[0].vars[3]->name=="Velocity_y"); - ASSERT(visData[0].vars[4]->name=="Velocity_z"); - Array& VelxData = visData[0].vars[2]->data; - Array& VelyData = visData[0].vars[3]->data; - Array& VelzData = visData[0].vars[4]->data; - fillData.copy(Averages.Vel_x,VelxData); - fillData.copy(Averages.Vel_y,VelyData); - fillData.copy(Averages.Vel_z,VelzData); - - ASSERT(visData[0].vars[6]->name=="BlobID"); - Array& BlobData = visData[0].vars[6]->data; - fillData.copy(Averages.Label_NWP,BlobData); + ASSERT( visData[0].vars[5]->name == "SignDist" ); + Array &SignData = visData[0].vars[5]->data; + fillData.copy( Averages.SDs, SignData ); + + ASSERT( visData[0].vars[1]->name == "Pressure" ); + Array &PressData = visData[0].vars[1]->data; + fillData.copy( Averages.Press, PressData ); + + ASSERT( visData[0].vars[2]->name == "Velocity_x" ); + ASSERT( visData[0].vars[3]->name == "Velocity_y" ); + ASSERT( visData[0].vars[4]->name == "Velocity_z" ); + Array &VelxData = visData[0].vars[2]->data; + Array &VelyData = visData[0].vars[3]->data; + Array &VelzData = visData[0].vars[4]->data; + fillData.copy( Averages.Vel_x, VelxData ); + fillData.copy( Averages.Vel_y, VelyData ); + fillData.copy( Averages.Vel_z, VelzData ); + + ASSERT( visData[0].vars[6]->name == "BlobID" ); + Array &BlobData = visData[0].vars[6]->data; + fillData.copy( Averages.Label_NWP, BlobData ); IO::writeData( timestep, visData, comm.comm ); - - PROFILE_STOP("Save Vis",1); + + PROFILE_STOP( "Save Vis", 1 ); }; + private: WriteVisWorkItem(); int timestep; - std::vector& visData; - TwoPhase& Averages; - fillHalo& fillData; + std::array n; + RankInfoStruct rank_info; + std::vector &visData; + TwoPhase &Averages; runAnalysis::commWrapper comm; }; // Helper class to write the vis file from a thread -class IOWorkItem: public ThreadPool::WorkItemRet +class IOWorkItem : public ThreadPool::WorkItemRet { public: - IOWorkItem(int timestep_, std::shared_ptr input_db_, std::vector& visData_, - SubPhase& Averages_, fillHalo& fillData_, runAnalysis::commWrapper&& comm_ ): - timestep(timestep_), input_db(input_db_), visData(visData_), Averages(Averages_), fillData(fillData_), comm(std::move(comm_)) - { - } - ~IOWorkItem() { } - virtual void run() { - auto color_db = input_db->getDatabase( "Color" ); - auto vis_db = input_db->getDatabase( "Visualization" ); - // int timestep = color_db->getWithDefault( "timestep", 0 ); + IOWorkItem( int timestep_, std::shared_ptr input_db_, + std::vector &visData_, SubPhase &Averages_, std::array n_, + RankInfoStruct rank_info_, runAnalysis::commWrapper &&comm_ ) + : timestep( timestep_ ), + input_db( input_db_ ), + visData( visData_ ), + Averages( Averages_ ), + n( std::move( n_ ) ), + rank_info( std::move( rank_info_ ) ), + comm( std::move( comm_ ) ) + { + } + ~IOWorkItem() {} + virtual void run() + { + PROFILE_START( "Save Vis", 1 ); - PROFILE_START("Save Vis",1); + auto color_db = input_db->getDatabase( "Color" ); + auto vis_db = input_db->getDatabase( "Visualization" ); + // int timestep = color_db->getWithDefault( "timestep", 0 ); - if (vis_db->getWithDefault( "save_phase_field", true )){ - ASSERT(visData[0].vars[0]->name=="phase"); - Array& PhaseData = visData[0].vars[0]->data; - fillData.copy(Averages.Phi,PhaseData); + fillHalo fillData( comm.comm, rank_info, n, { 1, 1, 1 }, 0, 1 ); + + if ( vis_db->getWithDefault( "save_phase_field", true ) ) { + ASSERT( visData[0].vars[0]->name == "phase" ); + Array &PhaseData = visData[0].vars[0]->data; + fillData.copy( Averages.Phi, PhaseData ); } - if (vis_db->getWithDefault( "save_pressure", false )){ - ASSERT(visData[0].vars[1]->name=="Pressure"); - Array& PressData = visData[0].vars[1]->data; - fillData.copy(Averages.Pressure,PressData); + if ( vis_db->getWithDefault( "save_pressure", false ) ) { + ASSERT( visData[0].vars[1]->name == "Pressure" ); + Array &PressData = visData[0].vars[1]->data; + fillData.copy( Averages.Pressure, PressData ); } - if (vis_db->getWithDefault( "save_velocity", false )){ - ASSERT(visData[0].vars[2]->name=="Velocity_x"); - ASSERT(visData[0].vars[3]->name=="Velocity_y"); - ASSERT(visData[0].vars[4]->name=="Velocity_z"); - Array& VelxData = visData[0].vars[2]->data; - Array& VelyData = visData[0].vars[3]->data; - Array& VelzData = visData[0].vars[4]->data; - fillData.copy(Averages.Vel_x,VelxData); - fillData.copy(Averages.Vel_y,VelyData); - fillData.copy(Averages.Vel_z,VelzData); + if ( vis_db->getWithDefault( "save_velocity", false ) ) { + ASSERT( visData[0].vars[2]->name == "Velocity_x" ); + ASSERT( visData[0].vars[3]->name == "Velocity_y" ); + ASSERT( visData[0].vars[4]->name == "Velocity_z" ); + Array &VelxData = visData[0].vars[2]->data; + Array &VelyData = visData[0].vars[3]->data; + Array &VelzData = visData[0].vars[4]->data; + fillData.copy( Averages.Vel_x, VelxData ); + fillData.copy( Averages.Vel_y, VelyData ); + fillData.copy( Averages.Vel_z, VelzData ); } - if (vis_db->getWithDefault( "save_distance", false )){ - ASSERT(visData[0].vars[5]->name=="SignDist"); - Array& SignData = visData[0].vars[5]->data; - fillData.copy(Averages.SDs,SignData); + if ( vis_db->getWithDefault( "save_distance", false ) ) { + ASSERT( visData[0].vars[5]->name == "SignDist" ); + Array &SignData = visData[0].vars[5]->data; + fillData.copy( Averages.SDs, SignData ); } - if (vis_db->getWithDefault( "save_connected_components", false )){ - ASSERT(visData[0].vars[6]->name=="BlobID"); - Array& BlobData = visData[0].vars[6]->data; - fillData.copy(Averages.morph_n->label,BlobData); - } - - if (vis_db->getWithDefault( "write_silo", true )) - IO::writeData( timestep, visData, comm.comm ); - - if (vis_db->getWithDefault( "save_8bit_raw", true )){ - char CurrentIDFilename[40]; - sprintf(CurrentIDFilename,"id_t%d.raw",timestep); - Averages.AggregateLabels(CurrentIDFilename); + if ( vis_db->getWithDefault( "save_connected_components", false ) ) { + ASSERT( visData[0].vars[6]->name == "BlobID" ); + Array &BlobData = visData[0].vars[6]->data; + fillData.copy( Averages.morph_n->label, BlobData ); } - PROFILE_STOP("Save Vis",1); + if ( vis_db->getWithDefault( "write_silo", true ) ) + IO::writeData( timestep, visData, comm.comm ); + + if ( vis_db->getWithDefault( "save_8bit_raw", true ) ) { + char CurrentIDFilename[40]; + sprintf( CurrentIDFilename, "id_t%d.raw", timestep ); + Averages.AggregateLabels( CurrentIDFilename ); + } + + PROFILE_STOP( "Save Vis", 1 ); }; + private: IOWorkItem(); int timestep; + std::array n; + RankInfoStruct rank_info; std::shared_ptr input_db; - std::vector& visData; - SubPhase& Averages; - fillHalo& fillData; + std::vector &visData; + SubPhase &Averages; runAnalysis::commWrapper comm; }; // Helper class to run the analysis from within a thread // Note: Averages will be modified after the constructor is called -class AnalysisWorkItem: public ThreadPool::WorkItemRet +class AnalysisWorkItem : public ThreadPool::WorkItemRet { public: - 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_) { } - ~AnalysisWorkItem() { } - virtual void run() { + 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_ ) + { + } + ~AnalysisWorkItem() {} + virtual void run() + { 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 ( matches(type,AnalysisType::CopyPhaseIndicator) ) { + Averages.Label_NWP = blob_ids->second; + Averages.Label_NWP_map = *id_list; + Averages.NumberComponents_WP = 1; + Averages.Label_WP.fill( 0.0 ); + if ( matches( type, AnalysisType::CopyPhaseIndicator ) ) { // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus); } - if ( matches(type,AnalysisType::ComputeAverages) ) { - PROFILE_START("Compute dist",1); + if ( matches( type, AnalysisType::ComputeAverages ) ) { + PROFILE_START( "Compute dist", 1 ); Averages.Initialize(); Averages.ComputeDelPhi(); - Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn); - Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus); - Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus); + Averages.ColorToSignedDistance( beta, Averages.Phase, Averages.SDn ); + Averages.ColorToSignedDistance( beta, Averages.Phase_tminus, Averages.Phase_tminus ); + Averages.ColorToSignedDistance( beta, Averages.Phase_tplus, Averages.Phase_tplus ); Averages.UpdateMeshValues(); Averages.ComputeLocal(); Averages.Reduce(); - Averages.PrintAll(timestep); + Averages.PrintAll( timestep ); Averages.Initialize(); Averages.ComponentAverages(); Averages.SortBlobs(); - Averages.PrintComponents(timestep); - PROFILE_STOP("Compute dist",1); + Averages.PrintComponents( timestep ); + PROFILE_STOP( "Compute dist", 1 ); } } + private: AnalysisWorkItem(); AnalysisType type; int timestep; - TwoPhase& Averages; + TwoPhase &Averages; BlobIDstruct blob_ids; BlobIDList id_list; double beta; }; -class TCATWorkItem: public ThreadPool::WorkItemRet +class TCATWorkItem : public ThreadPool::WorkItemRet { public: - TCATWorkItem( 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_) { } - ~TCATWorkItem() { } - virtual void run() { + TCATWorkItem( 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_ ) + { + } + ~TCATWorkItem() {} + virtual void run() + { 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 ( matches(type,AnalysisType::CopyPhaseIndicator) ) { + Averages.Label_NWP = blob_ids->second; + Averages.Label_NWP_map = *id_list; + Averages.NumberComponents_WP = 1; + Averages.Label_WP.fill( 0.0 ); + if ( matches( type, AnalysisType::CopyPhaseIndicator ) ) { // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus); } - if ( matches(type,AnalysisType::ComputeAverages) ) { - PROFILE_START("Compute TCAT",1); + if ( matches( type, AnalysisType::ComputeAverages ) ) { + PROFILE_START( "Compute TCAT", 1 ); Averages.Initialize(); Averages.ComputeDelPhi(); - Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn); - Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus); - Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus); + Averages.ColorToSignedDistance( beta, Averages.Phase, Averages.SDn ); + Averages.ColorToSignedDistance( beta, Averages.Phase_tminus, Averages.Phase_tminus ); + Averages.ColorToSignedDistance( beta, Averages.Phase_tplus, Averages.Phase_tplus ); Averages.UpdateMeshValues(); Averages.ComputeLocal(); Averages.Reduce(); - Averages.PrintAll(timestep); - PROFILE_STOP("Compute TCAT",1); + Averages.PrintAll( timestep ); + PROFILE_STOP( "Compute TCAT", 1 ); } } + private: TCATWorkItem(); AnalysisType type; int timestep; - TwoPhase& Averages; + TwoPhase &Averages; BlobIDstruct blob_ids; BlobIDList id_list; double beta; }; -class GanglionTrackingWorkItem: public ThreadPool::WorkItemRet +class GanglionTrackingWorkItem : public ThreadPool::WorkItemRet { public: - GanglionTrackingWorkItem( 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_) { } - ~GanglionTrackingWorkItem() { } - virtual void run() { + GanglionTrackingWorkItem( 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_ ) + { + } + ~GanglionTrackingWorkItem() {} + virtual void run() + { 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 ( matches(type,AnalysisType::CopyPhaseIndicator) ) { + Averages.Label_NWP = blob_ids->second; + Averages.Label_NWP_map = *id_list; + Averages.NumberComponents_WP = 1; + Averages.Label_WP.fill( 0.0 ); + if ( matches( type, AnalysisType::CopyPhaseIndicator ) ) { // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus); } - if ( matches(type,AnalysisType::ComputeAverages) ) { - PROFILE_START("Compute ganglion",1); + if ( matches( type, AnalysisType::ComputeAverages ) ) { + PROFILE_START( "Compute ganglion", 1 ); Averages.Initialize(); Averages.ComputeDelPhi(); - Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn); - Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus); - Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus); + Averages.ColorToSignedDistance( beta, Averages.Phase, Averages.SDn ); + Averages.ColorToSignedDistance( beta, Averages.Phase_tminus, Averages.Phase_tminus ); + Averages.ColorToSignedDistance( beta, Averages.Phase_tplus, Averages.Phase_tplus ); Averages.UpdateMeshValues(); Averages.ComponentAverages(); Averages.SortBlobs(); - Averages.PrintComponents(timestep); - PROFILE_STOP("Compute ganglion",1); + Averages.PrintComponents( timestep ); + PROFILE_STOP( "Compute ganglion", 1 ); } } + private: GanglionTrackingWorkItem(); AnalysisType type; int timestep; - TwoPhase& Averages; + TwoPhase &Averages; BlobIDstruct blob_ids; BlobIDList id_list; double beta; }; -class BasicWorkItem: public ThreadPool::WorkItemRet +class BasicWorkItem : public ThreadPool::WorkItemRet { public: - BasicWorkItem( AnalysisType type_, int timestep_, SubPhase& Averages_ ): - type(type_), timestep(timestep_), Averages(Averages_){ } - ~BasicWorkItem() { } - virtual void run() { + BasicWorkItem( AnalysisType type_, int timestep_, SubPhase &Averages_ ) + : type( type_ ), timestep( timestep_ ), Averages( Averages_ ) + { + } + ~BasicWorkItem() {} + virtual void run() + { - if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { + if ( matches( type, AnalysisType::CopyPhaseIndicator ) ) { // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus); } - if ( matches(type,AnalysisType::ComputeAverages) ) { - PROFILE_START("Compute basic averages",1); + if ( matches( type, AnalysisType::ComputeAverages ) ) { + PROFILE_START( "Compute basic averages", 1 ); Averages.Basic(); - PROFILE_STOP("Compute basic averages",1); + PROFILE_STOP( "Compute basic averages", 1 ); } } + private: BasicWorkItem(); AnalysisType type; int timestep; - SubPhase& Averages; + SubPhase &Averages; double beta; }; -class SubphaseWorkItem: public ThreadPool::WorkItemRet +class SubphaseWorkItem : public ThreadPool::WorkItemRet { public: - SubphaseWorkItem( AnalysisType type_, int timestep_, SubPhase& Averages_ ): - type(type_), timestep(timestep_), Averages(Averages_){ } - ~SubphaseWorkItem() { } - virtual void run() { - - PROFILE_START("Compute subphase",1); - Averages.Full(); - Averages.Write(timestep); - PROFILE_STOP("Compute subphase",1); + SubphaseWorkItem( AnalysisType type_, int timestep_, SubPhase &Averages_ ) + : type( type_ ), timestep( timestep_ ), Averages( Averages_ ) + { } + ~SubphaseWorkItem() {} + virtual void run() + { + + PROFILE_START( "Compute subphase", 1 ); + Averages.Full(); + Averages.Write( timestep ); + PROFILE_STOP( "Compute subphase", 1 ); + } + private: SubphaseWorkItem(); AnalysisType type; int timestep; - SubPhase& Averages; + SubPhase &Averages; double beta; }; - /****************************************************************** * MPI comm wrapper for use with analysis * ******************************************************************/ -runAnalysis::commWrapper::commWrapper( int tag_, const Utilities::MPI& comm_, runAnalysis* analysis_ ): - comm(comm_), - tag(tag_), - analysis(analysis_) +runAnalysis::commWrapper::commWrapper( + int tag_, const Utilities::MPI &comm_, runAnalysis *analysis_ ) + : comm( comm_ ), tag( tag_ ), analysis( analysis_ ) { } -runAnalysis::commWrapper::commWrapper( commWrapper &&rhs ): - comm(rhs.comm), - tag(rhs.tag), - analysis(rhs.analysis) +runAnalysis::commWrapper::commWrapper( commWrapper &&rhs ) + : comm( rhs.comm ), tag( rhs.tag ), analysis( rhs.analysis ) { rhs.tag = -1; } @@ -482,48 +563,44 @@ runAnalysis::commWrapper::~commWrapper() comm.barrier(); analysis->d_comm_used[tag] = false; } -runAnalysis::commWrapper runAnalysis::getComm( ) +runAnalysis::commWrapper runAnalysis::getComm() { // Get a tag from root int tag = -1; if ( d_rank == 0 ) { - for (int i=0; i<1024; i++) { + for ( int i = 0; i < 1024; i++ ) { if ( !d_comm_used[i] ) { tag = i; break; } } if ( tag == -1 ) - ERROR("Unable to get comm"); + ERROR( "Unable to get comm" ); } - tag = d_comm.bcast( tag, 0 ); + tag = d_comm.bcast( tag, 0 ); d_comm_used[tag] = true; if ( d_comms[tag].isNull() ) d_comms[tag] = d_comm.dup(); - return commWrapper(tag,d_comms[tag],this); + return commWrapper( tag, d_comms[tag], this ); } /****************************************************************** * Constructor/Destructors * ******************************************************************/ -runAnalysis::runAnalysis( std::shared_ptr input_db, - const RankInfoStruct& rank_info, - std::shared_ptr ScaLBL_Comm, - std::shared_ptr Dm, - int Np, - bool Regular, - IntArray Map ): - d_Np( Np ), - d_regular ( Regular), - d_rank_info( rank_info ), - d_Map( Map ), - d_comm( Dm->Comm.dup() ), - d_ScaLBL_Comm( ScaLBL_Comm) +runAnalysis::runAnalysis( std::shared_ptr input_db, const RankInfoStruct &rank_info, + std::shared_ptr ScaLBL_Comm, std::shared_ptr Dm, int Np, + bool Regular, IntArray Map ) + : d_Np( Np ), + d_regular( Regular ), + d_rank_info( rank_info ), + d_Map( Map ), + d_comm( Dm->Comm.dup() ), + d_ScaLBL_Comm( ScaLBL_Comm ) { - auto db = input_db->getDatabase( "Analysis" ); - auto vis_db = input_db->getDatabase( "Visualization" ); + auto db = input_db->getDatabase( "Analysis" ); + auto vis_db = input_db->getDatabase( "Visualization" ); // Ids of work items to use for dependencies ThreadPool::thread_id_t d_wait_blobID; @@ -533,117 +610,118 @@ runAnalysis::runAnalysis( std::shared_ptr input_db, ThreadPool::thread_id_t d_wait_subphase; char rankString[20]; - sprintf(rankString,"%05d",Dm->rank()); - d_n[0] = Dm->Nx-2; - d_n[1] = Dm->Ny-2; - d_n[2] = Dm->Nz-2; + sprintf( rankString, "%05d", Dm->rank() ); + d_n[0] = Dm->Nx - 2; + d_n[1] = Dm->Ny - 2; + d_n[2] = Dm->Nz - 2; d_N[0] = Dm->Nx; d_N[1] = Dm->Ny; d_N[2] = Dm->Nz; - - d_restart_interval = db->getScalar( "restart_interval" ); - d_analysis_interval = db->getScalar( "analysis_interval" ); + + d_restart_interval = db->getScalar( "restart_interval" ); + d_analysis_interval = db->getScalar( "analysis_interval" ); d_subphase_analysis_interval = INT_MAX; - d_visualization_interval = INT_MAX; - d_blobid_interval = INT_MAX; - if (db->keyExists( "blobid_interval" )){ - d_blobid_interval = db->getScalar( "blobid_interval" ); - } - if (db->keyExists( "visualization_interval" )){ - d_visualization_interval = db->getScalar( "visualization_interval" ); - } - if (db->keyExists( "subphase_analysis_interval" )){ - d_subphase_analysis_interval = db->getScalar( "subphase_analysis_interval" ); - } - + d_visualization_interval = INT_MAX; + d_blobid_interval = INT_MAX; + if ( db->keyExists( "blobid_interval" ) ) { + d_blobid_interval = db->getScalar( "blobid_interval" ); + } + if ( db->keyExists( "visualization_interval" ) ) { + d_visualization_interval = db->getScalar( "visualization_interval" ); + } + if ( db->keyExists( "subphase_analysis_interval" ) ) { + d_subphase_analysis_interval = db->getScalar( "subphase_analysis_interval" ); + } + auto restart_file = db->getScalar( "restart_file" ); - d_restartFile = restart_file + "." + rankString; - - + d_restartFile = restart_file + "." + rankString; + + d_rank = d_comm.getRank(); - writeIDMap(ID_map_struct(),0,id_map_filename); + writeIDMap( ID_map_struct(), 0, id_map_filename ); // Initialize IO for silo - IO::initialize("","silo","false"); - // Create the MeshDataStruct - d_meshData.resize(1); + IO::initialize( "", "silo", "false" ); + // Create the MeshDataStruct + d_meshData.resize( 1 ); d_meshData[0].meshName = "domain"; - d_meshData[0].mesh = std::make_shared( d_rank_info,d_n[0],d_n[1],d_n[2],Dm->Lx,Dm->Ly,Dm->Lz ); - auto PhaseVar = std::make_shared(); - auto PressVar = std::make_shared(); - auto VxVar = std::make_shared(); - auto VyVar = std::make_shared(); - auto VzVar = std::make_shared(); + d_meshData[0].mesh = std::make_shared( + d_rank_info, d_n[0], d_n[1], d_n[2], Dm->Lx, Dm->Ly, Dm->Lz ); + auto PhaseVar = std::make_shared(); + auto PressVar = std::make_shared(); + auto VxVar = std::make_shared(); + auto VyVar = std::make_shared(); + auto VzVar = std::make_shared(); auto SignDistVar = std::make_shared(); - auto BlobIDVar = std::make_shared(); - - if (vis_db->getWithDefault( "save_phase_field", true )){ + auto BlobIDVar = std::make_shared(); + + if ( vis_db->getWithDefault( "save_phase_field", true ) ) { PhaseVar->name = "phase"; PhaseVar->type = IO::VariableType::VolumeVariable; - PhaseVar->dim = 1; - PhaseVar->data.resize(d_n[0],d_n[1],d_n[2]); - d_meshData[0].vars.push_back(PhaseVar); + PhaseVar->dim = 1; + PhaseVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( PhaseVar ); } - if (vis_db->getWithDefault( "save_pressure", false )){ + if ( vis_db->getWithDefault( "save_pressure", false ) ) { PressVar->name = "Pressure"; PressVar->type = IO::VariableType::VolumeVariable; - PressVar->dim = 1; - PressVar->data.resize(d_n[0],d_n[1],d_n[2]); - d_meshData[0].vars.push_back(PressVar); + PressVar->dim = 1; + PressVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( PressVar ); } - if (vis_db->getWithDefault( "save_velocity", false )){ + if ( vis_db->getWithDefault( "save_velocity", false ) ) { VxVar->name = "Velocity_x"; VxVar->type = IO::VariableType::VolumeVariable; - VxVar->dim = 1; - VxVar->data.resize(d_n[0],d_n[1],d_n[2]); - d_meshData[0].vars.push_back(VxVar); + VxVar->dim = 1; + VxVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( VxVar ); VyVar->name = "Velocity_y"; VyVar->type = IO::VariableType::VolumeVariable; - VyVar->dim = 1; - VyVar->data.resize(d_n[0],d_n[1],d_n[2]); - d_meshData[0].vars.push_back(VyVar); + VyVar->dim = 1; + VyVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( VyVar ); VzVar->name = "Velocity_z"; VzVar->type = IO::VariableType::VolumeVariable; - VzVar->dim = 1; - VzVar->data.resize(d_n[0],d_n[1],d_n[2]); - d_meshData[0].vars.push_back(VzVar); + VzVar->dim = 1; + VzVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( VzVar ); } - if (vis_db->getWithDefault( "save_distance", false )){ + if ( vis_db->getWithDefault( "save_distance", false ) ) { SignDistVar->name = "SignDist"; SignDistVar->type = IO::VariableType::VolumeVariable; - SignDistVar->dim = 1; - SignDistVar->data.resize(d_n[0],d_n[1],d_n[2]); - d_meshData[0].vars.push_back(SignDistVar); + SignDistVar->dim = 1; + SignDistVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( SignDistVar ); } - if (vis_db->getWithDefault( "save_connected_components", false )){ + if ( vis_db->getWithDefault( "save_connected_components", false ) ) { BlobIDVar->name = "BlobID"; BlobIDVar->type = IO::VariableType::VolumeVariable; - BlobIDVar->dim = 1; - BlobIDVar->data.resize(d_n[0],d_n[1],d_n[2]); - d_meshData[0].vars.push_back(BlobIDVar); + BlobIDVar->dim = 1; + BlobIDVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( BlobIDVar ); } - + // Initialize the comms - for (int i=0; i<1024; i++) + for ( int i = 0; i < 1024; i++ ) d_comm_used[i] = false; // Initialize the threads int N_threads = db->getWithDefault( "N_threads", 4 ); - auto method = db->getWithDefault( "load_balance", "default" ); + auto method = db->getWithDefault( "load_balance", "default" ); createThreads( method, N_threads ); } -runAnalysis::~runAnalysis( ) +runAnalysis::~runAnalysis() { // Finish processing analysis finish(); } -void runAnalysis::finish( ) +void runAnalysis::finish() { - PROFILE_START("finish"); + PROFILE_START( "finish" ); // Wait for the work items to finish d_tpool.wait_pool_finished(); // Clear the wait ids @@ -654,23 +732,23 @@ void runAnalysis::finish( ) d_wait_restart.reset(); // Syncronize d_comm.barrier(); - PROFILE_STOP("finish"); + PROFILE_STOP( "finish" ); } /****************************************************************** * Set the thread affinities * ******************************************************************/ -void print( const std::vector& ids ) +void print( const std::vector &ids ) { if ( ids.empty() ) return; - printf("%i",ids[0]); - for (size_t i=1; i 0 ) - std::cerr << "Warning: Failed to start MPI with necessary thread support, errors may occur\n"; + std::cerr + << "Warning: Failed to start MPI with necessary thread support, errors may occur\n"; // Create the threads const auto cores = d_tpool.getProcessAffinity(); if ( N_threads == 0 ) { @@ -694,17 +773,17 @@ void runAnalysis::createThreads( const std::string& method, int N_threads ) int N = cores.size() - 1; d_tpool.setNumThreads( N ); d_tpool.setThreadAffinity( { cores[0] } ); - for ( int i=0; i input_db, TwoPhase& Averages, const double *Phi, - double *Pressure, double *Velocity, double *fq, double *Den) +void runAnalysis::run( int timestep, std::shared_ptr input_db, TwoPhase &Averages, + const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den ) { - int N = d_N[0]*d_N[1]*d_N[2]; + int N = d_N[0] * d_N[1] * d_N[2]; NULL_USE( N ); NULL_USE( Phi ); - - auto db = input_db->getDatabase( "Analysis" ); - //int timestep = db->getWithDefault( "timestep", 0 ); + + auto db = input_db->getDatabase( "Analysis" ); + // int timestep = db->getWithDefault( "timestep", 0 ); // Check which analysis steps we need to perform auto type = computeAnalysisType( timestep ); @@ -783,18 +861,18 @@ void runAnalysis::run(int timestep, std::shared_ptr input_db, TwoPhase finish(); } - PROFILE_START("run"); + PROFILE_START( "run" ); // Copy the appropriate variables to the host (so we can spawn new threads) ScaLBL_DeviceBarrier(); - PROFILE_START("Copy data to host",1); + PROFILE_START( "Copy data to host", 1 ); std::shared_ptr phase; /* if ( matches(type,AnalysisType::CopyPhaseIndicator) || matches(type,AnalysisType::ComputeAverages) || - matches(type,AnalysisType::CopySimState) || + matches(type,AnalysisType::CopySimState) || matches(type,AnalysisType::IdentifyBlobs) ) { - phase = std::shared_ptr(new DoubleArray(d_N[0],d_N[1],d_N[2])); + phase = std::make_shared(d_N[0],d_N[1],d_N[2]); //ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double)); // try 2 d_ScaLBL_Comm.RegulLayout(d_Map,Phi,Averages.Phase); // memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double)); @@ -820,133 +898,137 @@ void runAnalysis::run(int timestep, std::shared_ptr input_db, TwoPhase delete [] TmpDat; } */ - //if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { - if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) { - if (d_regular) - d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tplus); - else - ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double)); - //memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double)); - } - if ( timestep%d_analysis_interval == 0 ) { - if (d_regular) - d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tminus); - else - ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double)); - //memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double)); - } - //if ( matches(type,AnalysisType::CopySimState) ) { - if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) { - // Copy the members of Averages to the cpu (phase was copied above) - PROFILE_START("Copy-Pressure",1); - ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np); - //ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np); - ScaLBL_DeviceBarrier(); - PROFILE_STOP("Copy-Pressure",1); - PROFILE_START("Copy-Wait",1); - PROFILE_STOP("Copy-Wait",1); - PROFILE_START("Copy-State",1); - //memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double)); - if (d_regular) - d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase); + // if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { + if ( timestep % d_analysis_interval + 8 == d_analysis_interval ) { + if ( d_regular ) + d_ScaLBL_Comm->RegularLayout( d_Map, Phi, Averages.Phase_tplus ); else - ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double)); + ScaLBL_CopyToHost( Averages.Phase_tplus.data(), Phi, N * sizeof( double ) ); + // memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double)); + } + if ( timestep % d_analysis_interval == 0 ) { + if ( d_regular ) + d_ScaLBL_Comm->RegularLayout( d_Map, Phi, Averages.Phase_tminus ); + else + ScaLBL_CopyToHost( Averages.Phase_tminus.data(), Phi, N * sizeof( double ) ); + // memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double)); + } + // if ( matches(type,AnalysisType::CopySimState) ) { + if ( timestep % d_analysis_interval + 4 == d_analysis_interval ) { + // Copy the members of Averages to the cpu (phase was copied above) + PROFILE_START( "Copy-Pressure", 1 ); + ScaLBL_D3Q19_Pressure( fq, Pressure, d_Np ); + // ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np); + ScaLBL_DeviceBarrier(); + PROFILE_STOP( "Copy-Pressure", 1 ); + PROFILE_START( "Copy-Wait", 1 ); + PROFILE_STOP( "Copy-Wait", 1 ); + PROFILE_START( "Copy-State", 1 ); + // memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double)); + if ( d_regular ) + d_ScaLBL_Comm->RegularLayout( d_Map, Phi, Averages.Phase ); + else + ScaLBL_CopyToHost( Averages.Phase.data(), Phi, N * sizeof( double ) ); // copy other variables - d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press); - d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x); - d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y); - d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z); - PROFILE_STOP("Copy-State",1); + d_ScaLBL_Comm->RegularLayout( d_Map, Pressure, Averages.Press ); + d_ScaLBL_Comm->RegularLayout( d_Map, &Velocity[0], Averages.Vel_x ); + d_ScaLBL_Comm->RegularLayout( d_Map, &Velocity[d_Np], Averages.Vel_y ); + d_ScaLBL_Comm->RegularLayout( d_Map, &Velocity[2 * d_Np], Averages.Vel_z ); + PROFILE_STOP( "Copy-State", 1 ); } - std::shared_ptr cfq,cDen; - //if ( matches(type,AnalysisType::CreateRestart) ) { - if (timestep%d_restart_interval==0){ + std::shared_ptr cfq, cDen; + // if ( matches(type,AnalysisType::CreateRestart) ) { + if ( timestep % d_restart_interval == 0 ) { // Copy restart data to the CPU - cDen = std::shared_ptr(new double[2*d_Np],DeleteArray); - cfq = std::shared_ptr(new double[19*d_Np],DeleteArray); - ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double)); - ScaLBL_CopyToHost(cDen.get(),Den,2*d_Np*sizeof(double)); + cDen = make_shared_array( 2 * d_Np ); + cfq = make_shared_array( 19 * d_Np ); + ScaLBL_CopyToHost( cfq.get(), fq, 19 * d_Np * sizeof( double ) ); + ScaLBL_CopyToHost( cDen.get(), Den, 2 * d_Np * sizeof( double ) ); } - PROFILE_STOP("Copy data to host",1); + PROFILE_STOP( "Copy data to host", 1 ); // Spawn threads to do blob identification work - if ( matches(type,AnalysisType::IdentifyBlobs) ) { - phase = std::shared_ptr(new DoubleArray(d_N[0],d_N[1],d_N[2])); - if (d_regular) - d_ScaLBL_Comm->RegularLayout(d_Map,Phi,*phase); + if ( matches( type, AnalysisType::IdentifyBlobs ) ) { + phase = std::make_shared( d_N[0], d_N[1], d_N[2] ); + if ( d_regular ) + d_ScaLBL_Comm->RegularLayout( d_Map, Phi, *phase ); else - ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double)); + ScaLBL_CopyToHost( phase->data(), Phi, N * sizeof( double ) ); - BlobIDstruct new_index(new std::pair(0,IntArray())); - BlobIDstruct new_ids(new std::pair(0,IntArray())); - BlobIDList new_list(new std::vector()); - auto work1 = new BlobIdentificationWorkItem1(timestep,d_N[0],d_N[1],d_N[2],d_rank_info, - phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm()); - auto work2 = new BlobIdentificationWorkItem2(timestep,d_N[0],d_N[1],d_N[2],d_rank_info, - phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm()); - work1->add_dependency(d_wait_blobID); - work2->add_dependency(d_tpool.add_work(work1)); - d_wait_blobID = d_tpool.add_work(work2); - d_last_index = new_index; - d_last_ids = new_ids; + auto new_index = std::make_shared>( 0, IntArray() ); + auto new_ids = std::make_shared>( 0, IntArray() ); + auto new_list = std::make_shared>(); + auto work1 = new BlobIdentificationWorkItem1( timestep, d_N[0], d_N[1], d_N[2], d_rank_info, + phase, Averages.SDs, d_last_ids, new_index, new_ids, new_list, getComm() ); + auto work2 = new BlobIdentificationWorkItem2( timestep, d_N[0], d_N[1], d_N[2], d_rank_info, + phase, Averages.SDs, d_last_ids, new_index, new_ids, new_list, getComm() ); + work1->add_dependency( d_wait_blobID ); + work2->add_dependency( d_tpool.add_work( work1 ) ); + d_wait_blobID = d_tpool.add_work( work2 ); + d_last_index = new_index; + d_last_ids = new_ids; d_last_id_map = new_list; } // Spawn threads to do the analysis work - //if (timestep%d_restart_interval==0){ + // if (timestep%d_restart_interval==0){ // if ( matches(type,AnalysisType::ComputeAverages) ) { - if ( timestep%d_analysis_interval == 0 ) { - auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta); - work->add_dependency(d_wait_blobID); - work->add_dependency(d_wait_analysis); - work->add_dependency(d_wait_vis); // Make sure we are done using analysis before modifying - d_wait_analysis = d_tpool.add_work(work); + if ( timestep % d_analysis_interval == 0 ) { + auto work = + new AnalysisWorkItem( type, timestep, Averages, d_last_index, d_last_id_map, d_beta ); + work->add_dependency( d_wait_blobID ); + work->add_dependency( d_wait_analysis ); + work->add_dependency( d_wait_vis ); // Make sure we are done using analysis before modifying + d_wait_analysis = d_tpool.add_work( work ); } // Spawn a thread to write the restart file // if ( matches(type,AnalysisType::CreateRestart) ) { - if (timestep%d_restart_interval==0){ + if ( timestep % d_restart_interval == 0 ) { - if (d_rank==0) { - input_db->putScalar( "Restart", true ); - std::ofstream OutStream("Restart.db"); - input_db->print(OutStream, ""); - OutStream.close(); - } - // Write the restart file (using a seperate thread) - auto work = new WriteRestartWorkItem(d_restartFile.c_str(),cDen,cfq,d_Np); - work->add_dependency(d_wait_restart); - d_wait_restart = d_tpool.add_work(work); + if ( d_rank == 0 ) { + input_db->putScalar( "Restart", true ); + std::ofstream OutStream( "Restart.db" ); + input_db->print( OutStream, "" ); + OutStream.close(); + } + // Write the restart file (using a seperate thread) + auto work = new WriteRestartWorkItem( d_restartFile.c_str(), cDen, cfq, d_Np ); + work->add_dependency( d_wait_restart ); + d_wait_restart = d_tpool.add_work( work ); } // Save the results for visualization // if ( matches(type,AnalysisType::CreateRestart) ) { - if (timestep%d_restart_interval==0){ + if ( timestep % d_restart_interval == 0 ) { // Write the vis files - commWrapper comm = getComm(); - fillHalo fillData( comm.comm, d_rank_info, d_n, {1,1,1}, 0, 1 ); - auto work = new WriteVisWorkItem( timestep, d_meshData, Averages, fillData, std::move( comm ) ); - work->add_dependency(d_wait_blobID); - work->add_dependency(d_wait_analysis); - work->add_dependency(d_wait_vis); - d_wait_vis = d_tpool.add_work(work); + auto work = + new WriteVisWorkItem( timestep, d_meshData, Averages, d_n, d_rank_info, getComm() ); + work->add_dependency( d_wait_blobID ); + work->add_dependency( d_wait_analysis ); + work->add_dependency( d_wait_vis ); + d_wait_vis = d_tpool.add_work( work ); } - PROFILE_STOP("run"); + PROFILE_STOP( "run" ); } /****************************************************************** * Run the analysis * ******************************************************************/ -void runAnalysis::basic(int timestep, std::shared_ptr input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den) +void runAnalysis::basic( int timestep, std::shared_ptr input_db, SubPhase &Averages, + const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den ) { - int N = d_N[0]*d_N[1]*d_N[2]; + int Nx = d_N[0]; + int Ny = d_N[1]; + int Nz = d_N[2]; + int N = Nx * Ny * Nz; NULL_USE( N ); // Check which analysis steps we need to perform - auto color_db = input_db->getDatabase( "Color" ); - auto vis_db = input_db->getDatabase( "Visualization" ); + auto color_db = input_db->getDatabase( "Color" ); + auto vis_db = input_db->getDatabase( "Visualization" ); - //int timestep = color_db->getWithDefault( "timestep", 0 ); + // int timestep = color_db->getWithDefault( "timestep", 0 ); auto type = computeAnalysisType( timestep ); if ( type == AnalysisType::AnalyzeNone ) return; @@ -957,101 +1039,102 @@ void runAnalysis::basic(int timestep, std::shared_ptr input_db, SubPha finish(); } - PROFILE_START("basic"); + PROFILE_START( "basic" ); // Copy the appropriate variables to the host (so we can spawn new threads) ScaLBL_DeviceBarrier(); - PROFILE_START("Copy data to host",1); + PROFILE_START( "Copy data to host", 1 ); - //if ( matches(type,AnalysisType::CopySimState) ) { - if ( timestep%d_analysis_interval == 0 ) { + // if ( matches(type,AnalysisType::CopySimState) ) { + if ( timestep % d_analysis_interval == 0 ) { finish(); // can't copy if threads are still working on data // Copy the members of Averages to the cpu (phase was copied above) - PROFILE_START("Copy-Pressure",1); - ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np); - //ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np); + PROFILE_START( "Copy-Pressure", 1 ); + ScaLBL_D3Q19_Pressure( fq, Pressure, d_Np ); + // ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np); ScaLBL_DeviceBarrier(); - PROFILE_STOP("Copy-Pressure",1); - PROFILE_START("Copy-Wait",1); - PROFILE_STOP("Copy-Wait",1); - PROFILE_START("Copy-State",1); + PROFILE_STOP( "Copy-Pressure", 1 ); + PROFILE_START( "Copy-Wait", 1 ); + PROFILE_STOP( "Copy-Wait", 1 ); + PROFILE_START( "Copy-State", 1 ); /*if (d_regular) d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phi); else */ - ScaLBL_CopyToHost(Averages.Phi.data(),Phi,N*sizeof(double)); - // copy other variables - d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Pressure); - d_ScaLBL_Comm->RegularLayout(d_Map,&Den[0],Averages.Rho_n); - d_ScaLBL_Comm->RegularLayout(d_Map,&Den[d_Np],Averages.Rho_w); - d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x); - d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y); - d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z); - PROFILE_STOP("Copy-State",1); + ScaLBL_CopyToHost( Averages.Phi.data(), Phi, N * sizeof( double ) ); + // copy other variables + d_ScaLBL_Comm->RegularLayout( d_Map, Pressure, Averages.Pressure ); + d_ScaLBL_Comm->RegularLayout( d_Map, &Den[0], Averages.Rho_n ); + d_ScaLBL_Comm->RegularLayout( d_Map, &Den[d_Np], Averages.Rho_w ); + d_ScaLBL_Comm->RegularLayout( d_Map, &Velocity[0], Averages.Vel_x ); + d_ScaLBL_Comm->RegularLayout( d_Map, &Velocity[d_Np], Averages.Vel_y ); + d_ScaLBL_Comm->RegularLayout( d_Map, &Velocity[2 * d_Np], Averages.Vel_z ); + PROFILE_STOP( "Copy-State", 1 ); } - PROFILE_STOP("Copy data to host"); + PROFILE_STOP( "Copy data to host" ); // Spawn threads to do the analysis work - //if (timestep%d_restart_interval==0){ + // if (timestep%d_restart_interval==0){ // if ( matches(type,AnalysisType::ComputeAverages) ) { - if ( timestep%d_analysis_interval == 0 ) { - auto work = new BasicWorkItem(type,timestep,Averages); - work->add_dependency(d_wait_subphase); // Make sure we are done using analysis before modifying - work->add_dependency(d_wait_analysis); - work->add_dependency(d_wait_vis); - d_wait_analysis = d_tpool.add_work(work); - } - - if ( timestep%d_subphase_analysis_interval == 0 ) { - auto work = new SubphaseWorkItem(type,timestep,Averages); - work->add_dependency(d_wait_subphase); // Make sure we are done using analysis before modifying - work->add_dependency(d_wait_analysis); - work->add_dependency(d_wait_vis); - d_wait_subphase = d_tpool.add_work(work); + if ( timestep % d_analysis_interval == 0 ) { + auto work = new BasicWorkItem( type, timestep, Averages ); + work->add_dependency( + d_wait_subphase ); // Make sure we are done using analysis before modifying + work->add_dependency( d_wait_analysis ); + work->add_dependency( d_wait_vis ); + d_wait_analysis = d_tpool.add_work( work ); } - if (timestep%d_restart_interval==0){ - std::shared_ptr cfq,cDen; - // Copy restart data to the CPU - cDen = std::shared_ptr(new double[2*d_Np],DeleteArray); - cfq = std::shared_ptr(new double[19*d_Np],DeleteArray); - ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double)); - ScaLBL_CopyToHost(cDen.get(),Den,2*d_Np*sizeof(double)); - - if (d_rank==0) { - color_db->putScalar("timestep",timestep); - color_db->putScalar( "Restart", true ); - input_db->putDatabase("Color", color_db); - std::ofstream OutStream("Restart.db"); - input_db->print(OutStream, ""); - OutStream.close(); - - } - // Write the restart file (using a seperate thread) - auto work1 = new WriteRestartWorkItem(d_restartFile.c_str(),cDen,cfq,d_Np); - work1->add_dependency(d_wait_restart); - d_wait_restart = d_tpool.add_work(work1); - + if ( timestep % d_subphase_analysis_interval == 0 ) { + auto work = new SubphaseWorkItem( type, timestep, Averages ); + work->add_dependency( + d_wait_subphase ); // Make sure we are done using analysis before modifying + work->add_dependency( d_wait_analysis ); + work->add_dependency( d_wait_vis ); + d_wait_subphase = d_tpool.add_work( work ); } - - if (timestep%d_visualization_interval==0){ + + if ( timestep % d_restart_interval == 0 ) { + std::shared_ptr cfq, cDen; + // Copy restart data to the CPU + cDen = make_shared_array( 2 * d_Np ); + cfq = make_shared_array( 19 * d_Np ); + ScaLBL_CopyToHost( cfq.get(), fq, 19 * d_Np * sizeof( double ) ); + ScaLBL_CopyToHost( cDen.get(), Den, 2 * d_Np * sizeof( double ) ); + + if ( d_rank == 0 ) { + color_db->putScalar( "timestep", timestep ); + color_db->putScalar( "Restart", true ); + input_db->putDatabase( "Color", color_db ); + std::ofstream OutStream( "Restart.db" ); + input_db->print( OutStream, "" ); + OutStream.close(); + } + // Write the restart file (using a seperate thread) + auto work1 = new WriteRestartWorkItem( d_restartFile.c_str(), cDen, cfq, d_Np ); + work1->add_dependency( d_wait_restart ); + d_wait_restart = d_tpool.add_work( work1 ); + } + + if ( timestep % d_visualization_interval == 0 ) { // Write the vis files - commWrapper comm = getComm(); - fillHalo fillData( comm.comm, d_rank_info, d_n, {1,1,1}, 0, 1 ); - auto work = new IOWorkItem( timestep, input_db, d_meshData, Averages, fillData, std::move( comm ) ); - work->add_dependency(d_wait_analysis); - work->add_dependency(d_wait_subphase); - work->add_dependency(d_wait_vis); - d_wait_vis = d_tpool.add_work(work); + auto work = + new IOWorkItem( timestep, input_db, d_meshData, Averages, d_n, d_rank_info, getComm() ); + work->add_dependency( d_wait_analysis ); + work->add_dependency( d_wait_subphase ); + work->add_dependency( d_wait_vis ); + d_wait_vis = d_tpool.add_work( work ); } - PROFILE_STOP("basic"); + PROFILE_STOP( "basic" ); } -void runAnalysis::WriteVisData(int timestep, std::shared_ptr input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den) +void runAnalysis::WriteVisData( int timestep, std::shared_ptr input_db, + SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, + double *Den ) { - auto color_db = input_db->getDatabase( "Color" ); - auto vis_db = input_db->getDatabase( "Visualization" ); - //int timestep = color_db->getWithDefault( "timestep", 0 ); + auto color_db = input_db->getDatabase( "Color" ); + auto vis_db = input_db->getDatabase( "Visualization" ); + // int timestep = color_db->getWithDefault( "timestep", 0 ); // Check which analysis steps we need to perform auto type = computeAnalysisType( timestep ); @@ -1067,16 +1150,15 @@ void runAnalysis::WriteVisData(int timestep, std::shared_ptr input_db, // Copy the appropriate variables to the host (so we can spawn new threads) ScaLBL_DeviceBarrier(); - PROFILE_START("write vis",1); + PROFILE_START( "write vis", 1 ); // if (Averages.WriteVis == true){ - commWrapper comm = getComm(); - fillHalo fillData( comm.comm, d_rank_info, d_n, {1,1,1}, 0, 1 ); - auto work2 = new IOWorkItem(timestep, input_db, d_meshData, Averages, fillData, std::move( comm ) ); - work2->add_dependency(d_wait_vis); - d_wait_vis = d_tpool.add_work(work2); + auto work2 = + new IOWorkItem( timestep, input_db, d_meshData, Averages, d_n, d_rank_info, getComm() ); + work2->add_dependency( d_wait_vis ); + d_wait_vis = d_tpool.add_work( work2 ); - //Averages.WriteVis = false; - - PROFILE_STOP("write vis"); + // Averages.WriteVis = false; + + PROFILE_STOP( "write vis" ); } diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h index 33adbcb0..a82c4ba0 100644 --- a/analysis/runAnalysis.h +++ b/analysis/runAnalysis.h @@ -1,41 +1,48 @@ #ifndef RunAnalysis_H_INC #define RunAnalysis_H_INC -#include "analysis/analysis.h" -#include "analysis/TwoPhase.h" #include "analysis/SubPhase.h" +#include "analysis/TwoPhase.h" +#include "analysis/analysis.h" #include "common/Communication.h" #include "common/ScaLBL.h" #include "threadpool/thread_pool.h" #include -typedef std::shared_ptr> BlobIDstruct; -typedef std::shared_ptr> BlobIDList; - // Types of analysis -enum class AnalysisType : uint64_t { AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02, - CopySimState=0x04, ComputeAverages=0x08, CreateRestart=0x10, WriteVis=0x20, ComputeSubphase=0x40 }; +enum class AnalysisType : uint64_t { + AnalyzeNone = 0, + IdentifyBlobs = 0x01, + CopyPhaseIndicator = 0x02, + CopySimState = 0x04, + ComputeAverages = 0x08, + CreateRestart = 0x10, + WriteVis = 0x20, + ComputeSubphase = 0x40 +}; //! Class to run the analysis in multiple threads class runAnalysis { public: - //! Constructor - runAnalysis(std::shared_ptr db, const RankInfoStruct& rank_info, - std::shared_ptr ScaLBL_Comm, std::shared_ptr dm, int Np, bool Regular, IntArray Map ); + runAnalysis( std::shared_ptr db, const RankInfoStruct &rank_info, + std::shared_ptr ScaLBL_Comm, std::shared_ptr dm, int Np, + bool Regular, IntArray Map ); //! Destructor ~runAnalysis(); //! Run the next analysis - void run(int timestep, std::shared_ptr db, TwoPhase &Averages, const double *Phi, + void run( int timestep, std::shared_ptr db, TwoPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den ); - - void basic( int timestep, std::shared_ptr db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den ); - void WriteVisData(int timestep, std::shared_ptr vis_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den); + + void basic( int timestep, std::shared_ptr db, SubPhase &Averages, const double *Phi, + double *Pressure, double *Velocity, double *fq, double *Den ); + void WriteVisData( int timestep, std::shared_ptr vis_db, SubPhase &Averages, + const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den ); //! Finish all active analysis void finish(); @@ -44,7 +51,8 @@ public: * \brief Set the affinities * \details This function will create the analysis threads and set the affinity * of this thread and all analysis threads. If MPI_THREAD_MULTIPLE is not - * enabled, the analysis threads will be disabled and the analysis will run in the current thread. + * enabled, the analysis threads will be disabled and the analysis will run in the current + * thread. * @param[in] method Method used to control the affinities: * none - Don't use threads (runs all analysis in the current thread) * default - Create the specified number of threads, but don't load balance @@ -53,39 +61,36 @@ public: * that all threads run on independent cores * @param[in] N_threads Number of threads, only used by some of the methods */ - void createThreads( const std::string& method = "default", int N_threads = 4 ); + void createThreads( const std::string &method = "default", int N_threads = 4 ); private: - runAnalysis(); // Determine the analysis to perform AnalysisType computeAnalysisType( int timestep ); public: - class commWrapper { - public: + public: Utilities::MPI comm; int tag; runAnalysis *analysis; - commWrapper( int tag, const Utilities::MPI& comm, runAnalysis *analysis ); - commWrapper( ) = delete; + commWrapper( int tag, const Utilities::MPI &comm, runAnalysis *analysis ); + commWrapper() = delete; commWrapper( const commWrapper &rhs ) = delete; - commWrapper& operator=( const commWrapper &rhs ) = delete; + commWrapper &operator=( const commWrapper &rhs ) = delete; commWrapper( commWrapper &&rhs ); ~commWrapper(); }; // Get a comm (not thread safe) - commWrapper getComm( ); + commWrapper getComm(); private: - - std::array d_n; // Number of local cells - std::array d_N; // NNumber of local cells with ghosts + std::array d_n; // Number of local cells + std::array d_N; // Number of local cells with ghosts int d_Np; int d_rank; int d_restart_interval, d_analysis_interval, d_blobid_interval, d_visualization_interval; @@ -95,9 +100,9 @@ private: ThreadPool d_tpool; RankInfoStruct d_rank_info; IntArray d_Map; - BlobIDstruct d_last_ids; - BlobIDstruct d_last_index; - BlobIDList d_last_id_map; + std::shared_ptr> d_last_ids; + std::shared_ptr> d_last_index; + std::shared_ptr> d_last_id_map; std::vector d_meshData; std::string d_restartFile; Utilities::MPI d_comm; @@ -114,8 +119,6 @@ private: // Friends friend commWrapper::~commWrapper(); - }; #endif - diff --git a/common/Communication.h b/common/Communication.h index 4cd9ad70..5baaa962 100644 --- a/common/Communication.h +++ b/common/Communication.h @@ -67,6 +67,10 @@ public: //! Destructor ~fillHalo( ); + fillHalo() = delete; + fillHalo(const fillHalo&) = delete; + fillHalo& operator=(const fillHalo&) = delete; + /*! * @brief Communicate the halos * @param[in] array The array on which we fill the halos @@ -93,9 +97,6 @@ private: TYPE *mem; TYPE *send[3][3][3], *recv[3][3][3]; MPI_Request send_req[3][3][3], recv_req[3][3][3]; - fillHalo(); // Private empty constructor - fillHalo(const fillHalo&); // Private copy constructor - fillHalo& operator=(const fillHalo&); // Private assignment operator void pack( const Array& array, int i, int j, int k, TYPE *buffer ); void unpack( Array& array, int i, int j, int k, const TYPE *buffer ); }; diff --git a/common/Domain.cpp b/common/Domain.cpp index d552fb8a..ec24365d 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -558,17 +558,13 @@ void Domain::Decomp( const std::string& Filename ) int64_t z_transition_size = (nprocz*nz - (global_Nz - zStart))/2; if (z_transition_size < 0) z_transition_size=0; - char LocalRankFilename[40]; - char *loc_id; - loc_id = new char [(nx+2)*(ny+2)*(nz+2)]; - // Set up the sub-domains if (RANK==0){ printf("Distributing subdomains across %i processors \n",nprocs); printf("Process grid: %i x %i x %i \n",nprocx,nprocy,nprocz); printf("Subdomain size: %i x %i x %i \n",nx,ny,nz); printf("Size of transition region: %ld \n", z_transition_size); - + auto loc_id = new char [(nx+2)*(ny+2)*(nz+2)]; for (int kp=0; kpcommunicator = comm.communicator; this->comm_rank = comm.comm_rank; this->comm_size = comm.comm_size; - this->d_ranks = comm.d_ranks; this->d_isNull = comm.d_isNull; this->d_manage = comm.d_manage; this->d_maxTag = comm.d_maxTag; @@ -537,7 +530,6 @@ MPI_CLASS &MPI_CLASS::operator=( MPI_CLASS &&rhs ) std::swap( profile_level, rhs.profile_level ); std::swap( comm_rank, rhs.comm_rank ); std::swap( comm_size, rhs.comm_size ); - std::swap( d_ranks, rhs.d_ranks ); std::swap( d_maxTag, rhs.d_maxTag ); std::swap( d_currentTag, rhs.d_currentTag ); std::swap( d_count, rhs.d_count ); @@ -560,7 +552,6 @@ std::atomic_int d_global_count_self = { 1 }; MPI_CLASS::MPI_CLASS( MPI_Comm comm, bool manage ) { d_count = nullptr; - d_ranks = nullptr; d_manage = false; tmp_alignment = -1; // Check if we are using our version of comm_world @@ -623,11 +614,7 @@ MPI_CLASS::MPI_CLASS( MPI_Comm comm, bool manage ) } if ( d_manage ) ++N_MPI_Comm_created; - // Create d_ranks - if ( comm_size > 1 ) { - d_ranks = new int[comm_size]; - d_ranks[0] = -1; - } + #else // We are not using MPI, intialize based on the communicator NULL_USE( manage ); @@ -636,7 +623,7 @@ MPI_CLASS::MPI_CLASS( MPI_Comm comm, bool manage ) d_maxTag = mpi_max_tag; d_isNull = communicator == MPI_COMM_NULL; if ( d_isNull ) - comm_size = 0; + comm_size = 0; #endif if ( communicator == MPI_CLASS_COMM_WORLD ) { d_currentTag = d_global_currentTag_world1; @@ -663,34 +650,32 @@ MPI_CLASS::MPI_CLASS( MPI_Comm comm, bool manage ) ************************************************************************/ std::vector MPI_CLASS::globalRanks() const { - // Get my global rank if it has not been set - static int myGlobalRank = -1; - if ( myGlobalRank == -1 ) { -#ifdef USE_MPI - if ( MPI_active() ) - MPI_Comm_rank( MPI_CLASS_COMM_WORLD, &myGlobalRank ); -#else - myGlobalRank = 0; -#endif - } - // Check if we are dealing with a serial or null communicator - if ( comm_size == 1 ) - return std::vector( 1, myGlobalRank ); - if ( d_ranks == nullptr || communicator == MPI_COMM_NULL ) + if ( d_isNull ) return std::vector(); - // Fill d_ranks if necessary - if ( d_ranks[0] == -1 ) { - if ( communicator == MPI_CLASS_COMM_WORLD ) { - for ( int i = 0; i < comm_size; i++ ) - d_ranks[i] = i; - } else { - - MPI_ASSERT( myGlobalRank != -1 ); - this->allGather( myGlobalRank, d_ranks ); - } +#ifdef USE_MPI + // Get my global rank and size if it has not been set + static int globalRank = -1; + static int globalSize = -1; + if ( globalRank == -1 && MPI_active() ) { + MPI_Comm_rank( MPI_CLASS_COMM_WORLD, &globalRank ); + MPI_Comm_size( MPI_CLASS_COMM_WORLD, &globalSize ); } - // Return d_ranks - return std::vector( d_ranks, d_ranks + comm_size ); + // Check if we are dealing with a serial or global communicator + if ( comm_size == 1 ) + return std::vector( 1, globalRank ); + if ( comm_size == globalSize ) { + std::vector ranks( globalSize ); + for ( int i = 0; i < globalSize; i++ ) + ranks[i] = i; + return ranks; + } + // Get the global rank from each rank in the communicator + auto ranks = allGather( globalRank ); + std::sort( ranks.begin(), ranks.end() ); + return ranks; +#else + return std::vector( 1, 1 ); +#endif } @@ -2806,49 +2791,44 @@ MPI_Request MPI_CLASS::IrecvBytes( } - /************************************************************************ * sendrecv * ************************************************************************/ #if defined( USE_MPI ) template<> -void MPI_CLASS::sendrecv( const char* sendbuf, int sendcount, int dest, int sendtag, - char* recvbuf, int recvcount, int source, int recvtag ) const +void MPI_CLASS::sendrecv( const char *sendbuf, int sendcount, int dest, int sendtag, + char *recvbuf, int recvcount, int source, int recvtag ) const { PROFILE_START( "sendrecv", profile_level ); - MPI_Sendrecv( sendbuf, sendcount, MPI_CHAR, dest, sendtag, - recvbuf, recvcount, MPI_CHAR, source, recvtag, - communicator, MPI_STATUS_IGNORE ); + MPI_Sendrecv( sendbuf, sendcount, MPI_CHAR, dest, sendtag, recvbuf, recvcount, MPI_CHAR, source, + recvtag, communicator, MPI_STATUS_IGNORE ); PROFILE_STOP( "sendrecv", profile_level ); } template<> -void MPI_CLASS::sendrecv( const int* sendbuf, int sendcount, int dest, int sendtag, - int* recvbuf, int recvcount, int source, int recvtag ) const +void MPI_CLASS::sendrecv( const int *sendbuf, int sendcount, int dest, int sendtag, + int *recvbuf, int recvcount, int source, int recvtag ) const { PROFILE_START( "sendrecv", profile_level ); - MPI_Sendrecv( sendbuf, sendcount, MPI_INT, dest, sendtag, - recvbuf, recvcount, MPI_INT, source, recvtag, - communicator, MPI_STATUS_IGNORE ); + MPI_Sendrecv( sendbuf, sendcount, MPI_INT, dest, sendtag, recvbuf, recvcount, MPI_INT, source, + recvtag, communicator, MPI_STATUS_IGNORE ); PROFILE_STOP( "sendrecv", profile_level ); } template<> -void MPI_CLASS::sendrecv( const float* sendbuf, int sendcount, int dest, int sendtag, - float* recvbuf, int recvcount, int source, int recvtag ) const +void MPI_CLASS::sendrecv( const float *sendbuf, int sendcount, int dest, int sendtag, + float *recvbuf, int recvcount, int source, int recvtag ) const { PROFILE_START( "sendrecv", profile_level ); - MPI_Sendrecv( sendbuf, sendcount, MPI_FLOAT, dest, sendtag, - recvbuf, recvcount, MPI_FLOAT, source, recvtag, - communicator, MPI_STATUS_IGNORE ); + MPI_Sendrecv( sendbuf, sendcount, MPI_FLOAT, dest, sendtag, recvbuf, recvcount, MPI_FLOAT, + source, recvtag, communicator, MPI_STATUS_IGNORE ); PROFILE_STOP( "sendrecv", profile_level ); } template<> -void MPI_CLASS::sendrecv( const double* sendbuf, int sendcount, int dest, int sendtag, - double* recvbuf, int recvcount, int source, int recvtag ) const +void MPI_CLASS::sendrecv( const double *sendbuf, int sendcount, int dest, int sendtag, + double *recvbuf, int recvcount, int source, int recvtag ) const { PROFILE_START( "sendrecv", profile_level ); - MPI_Sendrecv( sendbuf, sendcount, MPI_DOUBLE, dest, sendtag, - recvbuf, recvcount, MPI_DOUBLE, source, recvtag, - communicator, MPI_STATUS_IGNORE ); + MPI_Sendrecv( sendbuf, sendcount, MPI_DOUBLE, dest, sendtag, recvbuf, recvcount, MPI_DOUBLE, + source, recvtag, communicator, MPI_STATUS_IGNORE ); PROFILE_STOP( "sendrecv", profile_level ); } #endif @@ -3815,17 +3795,16 @@ MPI MPI::loadBalance( double local, std::vector work ) MPI_ASSERT( (int) work.size() == getSize() ); auto perf = allGather( local ); std::vector I( work.size() ); - for ( size_t i=0; i key( work.size() ); - for ( size_t i=0; i globalRanks() const; @@ -796,7 +802,8 @@ public: // Member functions * @brief This function sends and recieves data using a blocking call */ template - void sendrecv( const type *sendbuf, int sendcount, int dest, int sendtag, type *recvbuf, int recvcount, int source, int recvtag ) const; + void sendrecv( const type *sendbuf, int sendcount, int dest, int sendtag, type *recvbuf, + int recvcount, int source, int recvtag ) const; /*! @@ -1126,9 +1133,6 @@ private: // data members // The rank and size of the communicator int comm_rank, comm_size; - // The ranks of the comm in the global comm - mutable int *volatile d_ranks; - // Some attributes int d_maxTag; int *volatile d_currentTag; diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index c3edc44f..182004ff 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -306,10 +306,130 @@ ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr Dm){ } -ScaLBL_Communicator::~ScaLBL_Communicator(){ - // destrutor does nothing (bad idea) - // -- note that there needs to be a way to free memory allocated on the device!!! +ScaLBL_Communicator::~ScaLBL_Communicator() +{ + + ScaLBL_FreeDeviceMemory( sendbuf_x ); + ScaLBL_FreeDeviceMemory( sendbuf_X ); + ScaLBL_FreeDeviceMemory( sendbuf_y ); + ScaLBL_FreeDeviceMemory( sendbuf_Y ); + ScaLBL_FreeDeviceMemory( sendbuf_z ); + ScaLBL_FreeDeviceMemory( sendbuf_Z ); + ScaLBL_FreeDeviceMemory( sendbuf_xy ); + ScaLBL_FreeDeviceMemory( sendbuf_xY ); + ScaLBL_FreeDeviceMemory( sendbuf_Xy ); + ScaLBL_FreeDeviceMemory( sendbuf_XY ); + ScaLBL_FreeDeviceMemory( sendbuf_xz ); + ScaLBL_FreeDeviceMemory( sendbuf_xZ ); + ScaLBL_FreeDeviceMemory( sendbuf_Xz ); + ScaLBL_FreeDeviceMemory( sendbuf_XZ ); + ScaLBL_FreeDeviceMemory( sendbuf_yz ); + ScaLBL_FreeDeviceMemory( sendbuf_yZ ); + ScaLBL_FreeDeviceMemory( sendbuf_Yz ); + ScaLBL_FreeDeviceMemory( sendbuf_YZ ); + ScaLBL_FreeDeviceMemory( recvbuf_x ); + ScaLBL_FreeDeviceMemory( recvbuf_X ); + ScaLBL_FreeDeviceMemory( recvbuf_y ); + ScaLBL_FreeDeviceMemory( recvbuf_Y ); + ScaLBL_FreeDeviceMemory( recvbuf_z ); + ScaLBL_FreeDeviceMemory( recvbuf_Z ); + ScaLBL_FreeDeviceMemory( recvbuf_xy ); + ScaLBL_FreeDeviceMemory( recvbuf_xY ); + ScaLBL_FreeDeviceMemory( recvbuf_Xy ); + ScaLBL_FreeDeviceMemory( recvbuf_XY ); + ScaLBL_FreeDeviceMemory( recvbuf_xz ); + ScaLBL_FreeDeviceMemory( recvbuf_xZ ); + ScaLBL_FreeDeviceMemory( recvbuf_Xz ); + ScaLBL_FreeDeviceMemory( recvbuf_XZ ); + ScaLBL_FreeDeviceMemory( recvbuf_yz ); + ScaLBL_FreeDeviceMemory( recvbuf_yZ ); + ScaLBL_FreeDeviceMemory( recvbuf_Yz ); + ScaLBL_FreeDeviceMemory( recvbuf_YZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_x ); + ScaLBL_FreeDeviceMemory( dvcSendList_X ); + ScaLBL_FreeDeviceMemory( dvcSendList_y ); + ScaLBL_FreeDeviceMemory( dvcSendList_Y ); + ScaLBL_FreeDeviceMemory( dvcSendList_z ); + ScaLBL_FreeDeviceMemory( dvcSendList_Z ); + ScaLBL_FreeDeviceMemory( dvcSendList_xy ); + ScaLBL_FreeDeviceMemory( dvcSendList_xY ); + ScaLBL_FreeDeviceMemory( dvcSendList_Xy ); + ScaLBL_FreeDeviceMemory( dvcSendList_XY ); + ScaLBL_FreeDeviceMemory( dvcSendList_xz ); + ScaLBL_FreeDeviceMemory( dvcSendList_xZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_Xz ); + ScaLBL_FreeDeviceMemory( dvcSendList_XZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_yz ); + ScaLBL_FreeDeviceMemory( dvcSendList_yZ ); + ScaLBL_FreeDeviceMemory( dvcSendList_Yz ); + ScaLBL_FreeDeviceMemory( dvcSendList_YZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_x ); + ScaLBL_FreeDeviceMemory( dvcRecvList_X ); + ScaLBL_FreeDeviceMemory( dvcRecvList_y ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Y ); + ScaLBL_FreeDeviceMemory( dvcRecvList_z ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Z ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xy ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xY ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Xy ); + ScaLBL_FreeDeviceMemory( dvcRecvList_XY ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_xZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Xz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_XZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_yz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_yZ ); + ScaLBL_FreeDeviceMemory( dvcRecvList_Yz ); + ScaLBL_FreeDeviceMemory( dvcRecvList_YZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_x ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_X ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_y ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Y ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_z ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Z ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xy ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xY ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Xy ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_XY ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_xZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Xz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_XZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_yz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_yZ ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_Yz ); + ScaLBL_FreeDeviceMemory( dvcRecvDist_YZ ); } +double ScaLBL_Communicator::GetPerformance(int *NeighborList, double *fq, int Np){ + /* EACH MPI PROCESS GETS ITS OWN MEASUREMENT*/ + /* use MRT kernels to check performance without communication / synchronization */ + int TIMESTEPS=500; + double RLX_SETA=1.0; + double RLX_SETB = 8.f*(2.f-RLX_SETA)/(8.f-RLX_SETA); + double FX = 0.0; + double FY = 0.0; + double FZ = 0.0; + ScaLBL_D3Q19_Init(fq, Np); + //.......create and start timer............ + double starttime,stoptime,cputime; + Barrier(); + starttime = MPI_Wtime(); + //......................................... + for (int t=0; t +#include +#include + +#define NBLOCKS 560 +#define NTHREADS 128 + +__global__ void dvc_ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradient, int start, int finish, int Np, int Nx, int Ny, int Nz) +{ + static int D3Q19[18][3]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}, + {1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0}, + {1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1}, + {0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}}; + + int i,j,k,n,N,idx; + int np,np2,nm; // neighbors + double v,vp,vp2,vm; // values at neighbors + double grad; + N = Nx*Ny*Nz; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s>>(Map, Phi, Gradient, start, finish, Np, Nx, Ny, Nz); + + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_MixedGradient: %s \n",cudaGetErrorString(err)); + } + cudaProfilerStop(); +} + diff --git a/example/Piston/input.db b/example/Piston/input.db index fab67cc5..5a9ba030 100644 --- a/example/Piston/input.db +++ b/example/Piston/input.db @@ -35,4 +35,6 @@ Analysis { load_balance = "independent" // Load balance method to use: "none", "default", "independent" } +Visualization { +} \ No newline at end of file diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 5f0e6f0e..3b8edd6c 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -10,14 +10,31 @@ color lattice boltzmann model #include ScaLBL_ColorModel::ScaLBL_ColorModel(int RANK, int NP, const Utilities::MPI& COMM): -rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rhoA(0),rhoB(0),alpha(0),beta(0), -Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0), -Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM) + rank(RANK), nprocs(NP), Restart(0), timestep(0), timestepMax(0), + tauA(0), tauB(0), rhoA(0), rhoB(0), alpha(0), beta(0), + Fx(0), Fy(0), Fz(0), flux(0), din(0), dout(0), + inletA(0), inletB(0), outletA(0), outletB(0), + Nx(0), Ny(0), Nz(0), N(0), Np(0), nprocx(0), nprocy(0), nprocz(0), + BoundaryCondition(0), Lx(0), Ly(0), Lz(0), id(nullptr), + NeighborList(nullptr), dvcMap(nullptr), fq(nullptr), Aq(nullptr), Bq(nullptr), + Den(nullptr), Phi(nullptr), ColorGrad(nullptr), Velocity(nullptr), Pressure(nullptr), + comm(COMM) { REVERSE_FLOW_DIRECTION = false; } -ScaLBL_ColorModel::~ScaLBL_ColorModel(){ - +ScaLBL_ColorModel::~ScaLBL_ColorModel() +{ + delete [] id; + ScaLBL_FreeDeviceMemory( NeighborList ); + ScaLBL_FreeDeviceMemory( dvcMap ); + ScaLBL_FreeDeviceMemory( fq ); + ScaLBL_FreeDeviceMemory( Aq ); + ScaLBL_FreeDeviceMemory( Bq ); + ScaLBL_FreeDeviceMemory( Den ); + ScaLBL_FreeDeviceMemory( Phi ); + ScaLBL_FreeDeviceMemory( Pressure ); + ScaLBL_FreeDeviceMemory( Velocity ); + ScaLBL_FreeDeviceMemory( ColorGrad ); } /*void ScaLBL_ColorModel::WriteCheckpoint(const char *FILENAME, const double *cPhi, const double *cfq, int Np) @@ -238,9 +255,26 @@ void ScaLBL_ColorModel::ReadInput(){ } } // MeanFilter(Averages->SDs); + Minkowski Solid(Dm); if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); CalcDist(Averages->SDs,id_solid,*Mask); - + Solid.ComputeScalar(Averages->SDs,0.0); + /* save averages */ + Averages->solid.V = Solid.Vi; + Averages->solid.A = Solid.Ai; + Averages->solid.H = Solid.Ji; + Averages->solid.X = Solid.Xi; + Averages->gsolid.V = Solid.Vi_global; + Averages->gsolid.A = Solid.Ai_global; + Averages->gsolid.H = Solid.Ji_global; + Averages->gsolid.X = Solid.Xi_global; + /* write to file */ + if (rank == 0) { + FILE *SOLID = fopen("solid.csv","w"); + fprintf(SOLID,"Vs As Hs Xs\n"); + fprintf(SOLID,"%.8g %.8g %.8g %.8g\n",Solid.Vi_global,Solid.Ai_global,Solid.Ji_global,Solid.Xi_global); + fclose(SOLID); + } if (rank == 0) cout << "Domain set." << endl; Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta); @@ -391,11 +425,13 @@ void ScaLBL_ColorModel::Create(){ // copy the neighbor list ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); + delete [] neighborList; // initialize phi based on PhaseLabel (include solid component labels) double *PhaseLabel; PhaseLabel = new double[N]; AssignComponentLabels(PhaseLabel); ScaLBL_CopyToDevice(Phi, PhaseLabel, N*sizeof(double)); + delete [] PhaseLabel; } /******************************************************** @@ -1080,7 +1116,6 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){ ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double)); // Extract only the connected part of NWP - BlobIDstruct new_index; double vF=0.0; double vS=0.0; ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm); comm.barrier(); @@ -1317,7 +1352,6 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta double volume_connected = 0.0; double second_biggest = 0.0; if (USE_CONNECTED_NWP){ - BlobIDstruct new_index; ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); comm.barrier(); diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index 3e19b717..01d13762 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -26,6 +26,8 @@ void ScaLBL_MRTModel::ReadParams(string filename){ tolerance = 1.0e-8; Fx = Fy = 0.0; Fz = 1.0e-5; + dout = 1.0; + din = 1.0; // Color Model parameters if (mrt_db->keyExists( "timestepMax" )){ @@ -194,7 +196,8 @@ void ScaLBL_MRTModel::Create(){ // copy the neighbor list ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); comm.barrier(); - + double MLUPS = ScaLBL_Comm->GetPerformance(NeighborList,fq,Np); + printf(" MLPUS=%f from rank %i\n",MLUPS,rank); } void ScaLBL_MRTModel::Initialize(){ diff --git a/sample_scripts/configure_huckleberry b/sample_scripts/configure_huckleberry index 8a4a313a..5ca6bb43 100755 --- a/sample_scripts/configure_huckleberry +++ b/sample_scripts/configure_huckleberry @@ -26,5 +26,5 @@ cmake \ -D USE_TIMER=0 \ ~/LBPM-WIA -make VERBOSE=1 -j8 && make install +make VERBOSE=1 -j4 && make install diff --git a/tests/TestNetcdf.cpp b/tests/TestNetcdf.cpp index 38fe08b3..6d43a04d 100644 --- a/tests/TestNetcdf.cpp +++ b/tests/TestNetcdf.cpp @@ -1,7 +1,7 @@ // Test reading/writing netcdf files #include "IO/netcdf.h" -#include "common/MPI_Helpers.h" +#include "common/MPI.h" #include "common/Communication.h" #include "common/UnitTest.h" @@ -13,7 +13,8 @@ void load( const std::string& ); void test_NETCDF( UnitTest& ut ) { - const int rank = comm_rank( MPI_COMM_WORLD ); + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); int nprocx = 2; int nprocy = 2; int nprocz = 2; @@ -26,11 +27,11 @@ void test_NETCDF( UnitTest& ut ) size_t z = info.kz*data.size(2); const char* filename = "test.nc"; std::vector dim = { (int) data.size(0)*nprocx, (int) data.size(1)*nprocy, (int) data.size(2)*nprocz }; - int fid = netcdf::open( filename, netcdf::CREATE, MPI_COMM_WORLD ); + int fid = netcdf::open( filename, netcdf::CREATE, comm ); auto dims = netcdf::defDim( fid, {"X", "Y", "Z"}, dim ); netcdf::write( fid, "tmp", dims, data, info ); netcdf::close( fid ); - MPI_Barrier( MPI_COMM_WORLD ); + comm.barrier(); // Read the contents of the file we created fid = netcdf::open( filename, netcdf::READ ); Array tmp = netcdf::getVar( fid, "tmp" ); diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index 6451d38a..590d5b8e 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -1,13 +1,13 @@ +#include +#include +#include +#include #include #include #include -#include -#include -#include -#include -#include "models/ColorModel.h" #include "common/Utilities.h" +#include "models/ColorModel.h" //#define WRE_SURFACES @@ -21,61 +21,59 @@ // Implementation of Two-Phase Immiscible LBM using CUDA //************************************************************************* -int main(int argc, char **argv) +int main( int argc, char **argv ) { - - // Initialize MPI - Utilities::startup( argc, argv ); - Utilities::MPI comm( MPI_COMM_WORLD ); - int rank = comm.getRank(); - int nprocs = comm.getSize(); - - // Load the input database - auto db = std::make_shared( argv[1] ); - // Initialize MPI and error handlers - auto multiple = db->getWithDefault( "MPI_THREAD_MULTIPLE", true ); - //Utilities::startup( argc, argv, multiple ); - //Utilities::MPI::changeProfileLevel( 1 ); + // Initialize + Utilities::startup( argc, argv ); - { // Limit scope so variables that contain communicators will free before MPI_Finialize + // Load the input database + auto db = std::make_shared( argv[1] ); - if (rank == 0){ - printf("********************************************************\n"); - printf("Running Color LBM \n"); - printf("********************************************************\n"); - } - // Initialize compute device - int device=ScaLBL_SetDevice(rank); - NULL_USE( device ); - ScaLBL_DeviceBarrier(); - comm.barrier(); + { // Limit scope so variables that contain communicators will free before MPI_Finialize - PROFILE_ENABLE(1); - //PROFILE_ENABLE_TRACE(); - //PROFILE_ENABLE_MEMORY(); - PROFILE_SYNCHRONIZE(); - PROFILE_START("Main"); - Utilities::setErrorHandlers(); + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); - auto filename = argv[1]; - ScaLBL_ColorModel ColorModel(rank,nprocs,comm); - ColorModel.ReadParams(filename); - ColorModel.SetDomain(); - ColorModel.ReadInput(); - ColorModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables - ColorModel.Initialize(); // initializing the model will set initial conditions for variables - ColorModel.Run(); - //ColorModel.WriteDebug(); + if ( rank == 0 ) { + printf( "********************************************************\n" ); + printf( "Running Color LBM \n" ); + printf( "********************************************************\n" ); + } + // Initialize compute device + int device = ScaLBL_SetDevice( rank ); + NULL_USE( device ); + ScaLBL_DeviceBarrier(); + comm.barrier(); - PROFILE_STOP("Main"); - auto file = db->getWithDefault( "TimerFile", "lbpm_color_simulator" ); - auto level = db->getWithDefault( "TimerLevel", 1 ); - PROFILE_SAVE(file,level); - // **************************************************** + PROFILE_ENABLE( 1 ); + // PROFILE_ENABLE_TRACE(); + // PROFILE_ENABLE_MEMORY(); + PROFILE_SYNCHRONIZE(); + PROFILE_START( "Main" ); + Utilities::setErrorHandlers(); + + auto filename = argv[1]; + ScaLBL_ColorModel ColorModel( rank, nprocs, comm ); + ColorModel.ReadParams( filename ); + ColorModel.SetDomain(); + ColorModel.ReadInput(); + ColorModel.Create(); // creating the model will create data structure to match the pore + // structure and allocate variables + ColorModel.Initialize(); // initializing the model will set initial conditions for variables + ColorModel.Run(); + // ColorModel.WriteDebug(); + + PROFILE_STOP( "Main" ); + auto file = db->getWithDefault( "TimerFile", "lbpm_color_simulator" ); + auto level = db->getWithDefault( "TimerLevel", 1 ); + PROFILE_SAVE( file, level ); + // **************************************************** - } // Limit scope so variables that contain communicators will free before MPI_Finialize + } // Limit scope so variables that contain communicators will free before MPI_Finialize - Utilities::shutdown(); + Utilities::shutdown(); + return 0; } diff --git a/tests/lbpm_morph_pp.cpp b/tests/lbpm_morph_pp.cpp index 12f6f319..e40dd6e0 100644 --- a/tests/lbpm_morph_pp.cpp +++ b/tests/lbpm_morph_pp.cpp @@ -128,7 +128,6 @@ int main(int argc, char **argv) comm.barrier(); // Extract only the connected part of NWP - BlobIDstruct new_index; double vF=0.0; double vS=0.0; ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,SignDist,vF,vS,phase_label,Dm->Comm); Dm->Comm.barrier(); diff --git a/tests/lbpm_uCT_pp.cpp b/tests/lbpm_uCT_pp.cpp index dbf9684b..b5d42e82 100644 --- a/tests/lbpm_uCT_pp.cpp +++ b/tests/lbpm_uCT_pp.cpp @@ -14,7 +14,7 @@ #include "common/Array.h" #include "common/Domain.h" #include "common/Communication.h" -#include "common/MPI_Helpers.h" +#include "common/MPI.h" #include "IO/MeshDatabase.h" #include "IO/Mesh.h" #include "IO/Writer.h" @@ -192,7 +192,7 @@ int main(int argc, char **argv) fillFloat[0]->fill( LOCVOL[0] ); } netcdf::close( fid ); - MPI_Barrier(comm); + comm.barrier(); PROFILE_STOP("ReadVolume"); if (rank==0) printf("Read complete\n"); @@ -255,15 +255,15 @@ int main(int argc, char **argv) } } } - count_plus=sumReduce( Dm[0]->Comm, count_plus); - count_minus=sumReduce( Dm[0]->Comm, count_minus); + count_plus = Dm[0]->Comm.sumReduce( count_plus); + count_minus = Dm[0]->Comm.sumReduce( count_minus); if (rank==0) printf("minimum value=%f, max value=%f \n",min_value,max_value); if (rank==0) printf("plus=%i, minus=%i \n",count_plus,count_minus); ASSERT( count_plus > 0 && count_minus > 0 ); - MPI_Barrier(comm); - mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / count_plus; - mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / count_minus; - MPI_Barrier(comm); + comm.barrier(); + mean_plus = Dm[0]->Comm.sumReduce( mean_plus ) / count_plus; + mean_minus = Dm[0]->Comm.sumReduce( mean_minus ) / count_minus; + comm.barrier(); if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus); //if (rank==0) printf("Scale the input data (size = %i) \n",LOCVOL[0].length()); @@ -284,7 +284,7 @@ int main(int argc, char **argv) // Fill the source data for the coarse meshes if (rank==0) printf("Coarsen the mesh for N_levels=%i \n",N_levels); - MPI_Barrier(comm); + comm.barrier(); PROFILE_START("CoarsenMesh"); for (int i=1; i filter(ratio[0],ratio[1],ratio[2]); @@ -300,7 +300,7 @@ int main(int argc, char **argv) printf(" filter_x=%i, filter_y=%i, filter_z=%i \n",int(filter.size(0)),int(filter.size(1)),int(filter.size(2)) ); printf(" ratio= %i,%i,%i \n",int(ratio[0]),int(ratio[1]),int(ratio[2]) ); } - MPI_Barrier(comm); + comm.barrier(); } PROFILE_STOP("CoarsenMesh"); @@ -312,7 +312,7 @@ int main(int argc, char **argv) NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx, rough_cutoff, lamda, nlm_sigsq, nlm_depth); PROFILE_STOP("Solve coarse mesh"); - MPI_Barrier(comm); + comm.barrier(); // Refine the solution PROFILE_START("Refine distance"); @@ -326,7 +326,7 @@ int main(int argc, char **argv) rough_cutoff, lamda, nlm_sigsq, nlm_depth); } PROFILE_STOP("Refine distance"); - MPI_Barrier(comm); + comm.barrier(); // Perform a final filter PROFILE_START("Filtering final domains"); @@ -424,14 +424,14 @@ int main(int argc, char **argv) meshData[0].vars.push_back(filter_Dist2_var); fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data ); #endif - MPI_Barrier(comm); + comm.barrier(); if (rank==0) printf("Writing output \n"); // Write visulization data IO::writeData( 0, meshData, comm ); if (rank==0) printf("Finished. \n"); // Compute the Minkowski functionals - MPI_Barrier(comm); + comm.barrier(); auto Averages = std::make_shared(Dm[0]); Array phase_label(Nx[0]+2,Ny[0]+2,Nz[0]+2);