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

Fixing bug with analysis

Browse Source
This commit is contained in:
Mark Berrill 2021-01-22 10:05:16 -05:00
parent c350e2b45e
commit 90ba7ed65a
15 changed files with 1229 additions and 926 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
.clang-format Normal file
View File

@ -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
...

375
IO/silo.hpp
View File

@ -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 <silo.h>
namespace silo {
/****************************************************
* Helper functions *
****************************************************/
template<class TYPE> static constexpr int getType();
template<> constexpr int getType<double>() { return DB_DOUBLE; }
template<> constexpr int getType<float>() { return DB_FLOAT; }
template<> constexpr int getType<int>() { return DB_INT; }
* Helper functions *
****************************************************/
template<class TYPE>
inline void copyData( Array<TYPE>& data, int type, const void *src )
static constexpr int getType();
template<>
constexpr int getType<double>()
{
return DB_DOUBLE;
}
template<>
constexpr int getType<float>()
{
return DB_FLOAT;
}
template<>
constexpr int getType<int>()
{
return DB_INT;
}
template<class TYPE>
inline void copyData( Array<TYPE> &data, int type, const void *src )
{
if ( type == getType<TYPE>() )
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<const double*>(src) );
data.copy( static_cast<const double *>( src ) );
else if ( type == DB_FLOAT )
data.copy( static_cast<const float*>(src) );
data.copy( static_cast<const float *>( src ) );
else if ( type == DB_INT )
data.copy( static_cast<const int*>(src) );
data.copy( static_cast<const int *>( src ) );
else
ERROR("Unknown type");
ERROR( "Unknown type" );
}
/****************************************************
* Write/read an arbitrary vector *
****************************************************/
template<class TYPE> constexpr int getSiloType();
template<> constexpr int getSiloType<int>() { return DB_INT; }
template<> constexpr int getSiloType<float>() { return DB_FLOAT; }
template<> constexpr int getSiloType<double>() { return DB_DOUBLE; }
* Write/read an arbitrary vector *
****************************************************/
template<class TYPE>
void write( DBfile* fid, const std::string& varname, const std::vector<TYPE>& data )
constexpr int getSiloType();
template<>
constexpr int getSiloType<int>()
{
return DB_INT;
}
template<>
constexpr int getSiloType<float>()
{
return DB_FLOAT;
}
template<>
constexpr int getSiloType<double>()
{
return DB_DOUBLE;
}
template<class TYPE>
void write( DBfile *fid, const std::string &varname, const std::vector<TYPE> &data )
{
int dims = data.size();
int err = DBWrite( fid, varname.c_str(), (void*) data.data(), &dims, 1, getSiloType<TYPE>() );
int err = DBWrite( fid, varname.c_str(), (void *) data.data(), &dims, 1, getSiloType<TYPE>() );
ASSERT( err == 0 );
}
template<class TYPE>
std::vector<TYPE> read( DBfile* fid, const std::string& varname )
std::vector<TYPE> read( DBfile *fid, const std::string &varname )
{
int N = DBGetVarLength( fid, varname.c_str() );
std::vector<TYPE> data(N);
std::vector<TYPE> data( N );
int err = DBReadVar( fid, varname.c_str(), data.data() );
ASSERT( err == 0 );
return data;
@ -66,31 +91,31 @@ std::vector<TYPE> read( DBfile* fid, const std::string& varname )
/****************************************************
* Helper function to get variable suffixes *
****************************************************/
* Helper function to get variable suffixes *
****************************************************/
inline std::vector<std::string> getVarSuffix( int ndim, int nvars )
{
std::vector<std::string> suffix(nvars);
std::vector<std::string> 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<std::string> getVarSuffix( int ndim, int nvars )
suffix[7] = "_zy";
suffix[8] = "_zz";
} else {
ERROR("Not finished");
ERROR( "Not finished" );
}
} else {
for (int i=0; i<nvars; i++)
suffix[i] = "_" + std::to_string(i+1);
for ( int i = 0; i < nvars; i++ )
suffix[i] = "_" + std::to_string( i + 1 );
}
return suffix;
}
/****************************************************
* Write/read a uniform mesh to silo *
****************************************************/
* Write/read a uniform mesh to silo *
****************************************************/
template<int NDIM>
void writeUniformMesh( DBfile* fid, const std::string& meshname,
const std::array<double,2*NDIM>& range, const std::array<int,NDIM>& N )
void writeUniformMesh( DBfile *fid, const std::string &meshname,
const std::array<double, 2 * NDIM> &range, const std::array<int, NDIM> &N )
{
PROFILE_START("writeUniformMesh",2);
PROFILE_START( "writeUniformMesh", 2 );
int dims[NDIM];
for (size_t d=0; d<N.size(); d++)
dims[d] = N[d]+1;
for ( size_t d = 0; d < N.size(); d++ )
dims[d] = N[d] + 1;
float *x = nullptr;
if ( NDIM >= 1 ) {
x = new float[dims[0]];
for (int i=0; i<N[0]; i++)
x[i] = range[0] + i*(range[1]-range[0])/N[0];
for ( int i = 0; i < N[0]; i++ )
x[i] = range[0] + i * ( range[1] - range[0] ) / N[0];
x[N[0]] = range[1];
}
float *y = nullptr;
if ( NDIM >= 2 ) {
y = new float[dims[1]];
for (int i=0; i<N[1]; i++)
y[i] = range[2] + i*(range[3]-range[2])/N[1];
for ( int i = 0; i < N[1]; i++ )
y[i] = range[2] + i * ( range[3] - range[2] ) / N[1];
y[N[1]] = range[3];
}
float *z = nullptr;
if ( NDIM >= 3 ) {
z = new float[dims[2]];
for (int i=0; i<N[2]; i++)
z[i] = range[4] + i*(range[5]-range[4])/N[2];
for ( int i = 0; i < N[2]; i++ )
z[i] = range[4] + i * ( range[5] - range[4] ) / N[2];
z[N[2]] = range[5];
}
float *coords[] = { x, y, z };
int err = DBPutQuadmesh( fid, meshname.c_str(), nullptr, coords, dims, NDIM, DB_FLOAT, DB_COLLINEAR, nullptr );
int err = DBPutQuadmesh(
fid, meshname.c_str(), nullptr, coords, dims, NDIM, DB_FLOAT, DB_COLLINEAR, nullptr );
delete[] x;
delete[] y;
delete[] z;
ASSERT( err == 0 );
PROFILE_STOP("writeUniformMesh",2);
PROFILE_STOP( "writeUniformMesh", 2 );
}
/****************************************************
* Write a vector/tensor quad variable *
****************************************************/
template<int NDIM,class TYPE>
void writeUniformMeshVariable( DBfile* fid, const std::string& meshname, const std::array<int,NDIM>& N,
const std::string& varname, const Array<TYPE>& data, VariableType type )
* Write a vector/tensor quad variable *
****************************************************/
template<int NDIM, class TYPE>
void writeUniformMeshVariable( DBfile *fid, const std::string &meshname,
const std::array<int, NDIM> &N, const std::string &varname, const Array<TYPE> &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<NDIM; d++)
ASSERT(N[d]+1==(int)data.size(d));
vartype = DB_NODECENT;
nvars = data.size(NDIM);
size_t N = data.length()/nvars;
for (int d=0; d<NDIM; d++)
dims[d] = data.size(d);
for (int i=0; i<nvars; i++)
vars[i] = &data(i*N);
ASSERT( data.ndim() == NDIM || data.ndim() == NDIM + 1 );
for ( int d = 0; d < NDIM; d++ )
ASSERT( N[d] + 1 == (int) data.size( d ) );
vartype = DB_NODECENT;
nvars = data.size( NDIM );
size_t N = data.length() / nvars;
for ( int d = 0; d < NDIM; d++ )
dims[d] = data.size( d );
for ( int i = 0; i < nvars; i++ )
vars[i] = &data( i * N );
} else if ( type == VariableType::EdgeVariable ) {
ERROR("Not finished");
ERROR( "Not finished" );
} else if ( type == VariableType::SurfaceVariable ) {
ERROR("Not finished");
ERROR( "Not finished" );
} else if ( type == VariableType::VolumeVariable ) {
ASSERT( data.ndim()==NDIM || data.ndim()==NDIM+1 );
for (int d=0; d<NDIM; d++)
ASSERT(N[d]==(int)data.size(d));
vartype = DB_ZONECENT;
nvars = data.size(NDIM);
size_t N = data.length()/nvars;
for (int d=0; d<NDIM; d++)
dims[d] = data.size(d);
for (int i=0; i<nvars; i++)
vars[i] = &data(i*N);
ASSERT( data.ndim() == NDIM || data.ndim() == NDIM + 1 );
for ( int d = 0; d < NDIM; d++ )
ASSERT( N[d] == (int) data.size( d ) );
vartype = DB_ZONECENT;
nvars = data.size( NDIM );
size_t N = data.length() / nvars;
for ( int d = 0; d < NDIM; d++ )
dims[d] = data.size( d );
for ( int i = 0; i < nvars; i++ )
vars[i] = &data( i * N );
} else {
ERROR("Invalid variable type");
ERROR( "Invalid variable type" );
}
auto suffix = getVarSuffix( NDIM, nvars );
std::vector<std::string> var_names(nvars);
for (int i=0; i<nvars; i++)
std::vector<std::string> var_names( nvars );
for ( int i = 0; i < nvars; i++ )
var_names[i] = varname + suffix[i];
std::vector<char*> varnames(nvars,nullptr);
for (int i=0; i<nvars; i++)
varnames[i] = const_cast<char*>(var_names[i].c_str());
int err = DBPutQuadvar( fid, varname.c_str(), meshname.c_str(), nvars,
varnames.data(), vars, dims, NDIM, nullptr, 0, getType<TYPE>(), vartype, nullptr );
std::vector<char *> varnames( nvars, nullptr );
for ( int i = 0; i < nvars; i++ )
varnames[i] = const_cast<char *>( var_names[i].c_str() );
int err = DBPutQuadvar( fid, varname.c_str(), meshname.c_str(), nvars, varnames.data(), vars,
dims, NDIM, nullptr, 0, getType<TYPE>(), vartype, nullptr );
ASSERT( err == 0 );
PROFILE_STOP("writeUniformMeshVariable",2);
PROFILE_STOP( "writeUniformMeshVariable", 2 );
}
template<class TYPE>
Array<TYPE> readUniformMeshVariable( DBfile* fid, const std::string& varname )
template<class TYPE>
Array<TYPE> readUniformMeshVariable( DBfile *fid, const std::string &varname )
{
auto var = DBGetQuadvar( fid, varname.c_str() );
ASSERT( var != nullptr );
Array<TYPE> data( var->nels, var->nvals );
int type = var->datatype;
for (int i=0; i<var->nvals; i++) {
for ( int i = 0; i < var->nvals; i++ ) {
Array<TYPE> data2( var->nels );
copyData<TYPE>( 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<size_t> dims( var->ndims+1, var->nvals );
for (int d=0; d<var->ndims; d++)
std::vector<size_t> 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<TYPE> readUniformMeshVariable( DBfile* fid, const std::string& varname )
/****************************************************
* Read/write a point mesh/variable to silo *
****************************************************/
* Read/write a point mesh/variable to silo *
****************************************************/
template<class TYPE>
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<TYPE>(), nullptr );
ASSERT( err == 0 );
}
template<class TYPE>
Array<TYPE> readPointMesh( DBfile* fid, const std::string& meshname )
template<class TYPE>
Array<TYPE> readPointMesh( DBfile *fid, const std::string &meshname )
{
auto mesh = DBGetPointmesh( fid, meshname.c_str() );
int N = mesh->nels;
int ndim = mesh->ndims;
Array<TYPE> coords(N,ndim);
int N = mesh->nels;
int ndim = mesh->ndims;
Array<TYPE> coords( N, ndim );
int type = mesh->datatype;
for (int d=0; d<ndim; d++) {
for ( int d = 0; d < ndim; d++ ) {
Array<TYPE> data2( N );
copyData<TYPE>( 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<class TYPE>
void writePointMeshVariable( DBfile* fid, const std::string& meshname,
const std::string& varname, const Array<TYPE>& data )
void writePointMeshVariable(
DBfile *fid, const std::string &meshname, const std::string &varname, const Array<TYPE> &data )
{
int N = data.size(0);
int nvars = data.size(1);
std::vector<const TYPE*> 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<TYPE>(), nullptr );
int N = data.size( 0 );
int nvars = data.size( 1 );
std::vector<const TYPE *> 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<TYPE>(), nullptr );
ASSERT( err == 0 );
}
template<class TYPE>
Array<TYPE> readPointMeshVariable( DBfile* fid, const std::string& varname )
template<class TYPE>
Array<TYPE> readPointMeshVariable( DBfile *fid, const std::string &varname )
{
auto var = DBGetPointvar( fid, varname.c_str() );
ASSERT( var != nullptr );
Array<TYPE> data( var->nels, var->nvals );
int type = var->datatype;
for (int i=0; i<var->nvals; i++) {
for ( int i = 0; i < var->nvals; i++ ) {
Array<TYPE> data2( var->nels );
copyData<TYPE>( 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<TYPE> readPointMeshVariable( DBfile* fid, const std::string& varname )
/****************************************************
* Read/write a triangle mesh *
****************************************************/
* Read/write a triangle mesh *
****************************************************/
template<class TYPE>
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<int> 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++)
std::vector<int> 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<TYPE>(), 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<TYPE>(), nullptr );
}
template<class TYPE>
void readTriMesh( DBfile* fid, const std::string& meshname, Array<TYPE>& coords, Array<int>& tri )
void readTriMesh( DBfile *fid, const std::string &meshname, Array<TYPE> &coords, Array<int> &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<ndim; d++) {
for ( int d = 0; d < ndim; d++ ) {
Array<TYPE> data2( N_nodes );
copyData<TYPE>( 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; i<N_zones; i++) {
for (int j=0; j<shapesize; j++)
tri(i,j) = zones->nodelist[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<class TYPE>
void writeTriMeshVariable( DBfile* fid, int ndim, const std::string& meshname,
const std::string& varname, const Array<TYPE>& data, VariableType type )
void writeTriMeshVariable( DBfile *fid, int ndim, const std::string &meshname,
const std::string &varname, const Array<TYPE> &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<nvars; i++)
vars[i] = &data(0,i);
nvars = data.size( 1 );
for ( int i = 0; i < nvars; i++ )
vars[i] = &data( 0, i );
} else if ( type == VariableType::EdgeVariable ) {
ERROR("Not finished");
ERROR( "Not finished" );
} else if ( type == VariableType::SurfaceVariable ) {
ERROR("Not finished");
ERROR( "Not finished" );
} else if ( type == VariableType::VolumeVariable ) {
vartype = DB_ZONECENT;
nvars = data.size(1);
for (int i=0; i<nvars; i++)
vars[i] = &data(0,i);
nvars = data.size( 1 );
for ( int i = 0; i < nvars; i++ )
vars[i] = &data( 0, i );
} else {
ERROR("Invalid variable type");
ERROR( "Invalid variable type" );
}
auto suffix = getVarSuffix( ndim, nvars );
std::vector<std::string> var_names(nvars);
for (int i=0; i<nvars; i++)
std::vector<std::string> var_names( nvars );
for ( int i = 0; i < nvars; i++ )
var_names[i] = varname + suffix[i];
std::vector<char*> varnames(nvars,nullptr);
for (int i=0; i<nvars; i++)
varnames[i] = const_cast<char*>(var_names[i].c_str());
DBPutUcdvar( fid, varname.c_str(), meshname.c_str(), nvars,
varnames.data(), vars, data.size(0), nullptr, 0, getType<TYPE>(), vartype, nullptr );
std::vector<char *> varnames( nvars, nullptr );
for ( int i = 0; i < nvars; i++ )
varnames[i] = const_cast<char *>( var_names[i].c_str() );
DBPutUcdvar( fid, varname.c_str(), meshname.c_str(), nvars, varnames.data(), vars,
data.size( 0 ), nullptr, 0, getType<TYPE>(), vartype, nullptr );
}
template<class TYPE>
Array<TYPE> readTriMeshVariable( DBfile* fid, const std::string& varname )
Array<TYPE> readTriMeshVariable( DBfile *fid, const std::string &varname )
{
auto var = DBGetUcdvar( fid, varname.c_str() );
ASSERT( var != nullptr );
Array<TYPE> data( var->nels, var->nvals );
int type = var->datatype;
for (int i=0; i<var->nvals; i++) {
for ( int i = 0; i < var->nvals; i++ ) {
Array<TYPE> data2( var->nels );
copyData<TYPE>( 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

1
analysis/morphology.cpp
View File

@ -542,7 +542,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
}
// 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();

1175
analysis/runAnalysis.cpp
View File

File diff suppressed because it is too large Load Diff

65
analysis/runAnalysis.h
View File

@ -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 <limits.h>
typedef std::shared_ptr<std::pair<int,IntArray>> BlobIDstruct;
typedef std::shared_ptr<std::vector<BlobIDType>> 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<Database> db, const RankInfoStruct& rank_info,
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> dm, int Np, bool Regular, IntArray Map );
runAnalysis( std::shared_ptr<Database> db, const RankInfoStruct &rank_info,
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr<Domain> dm, int Np,
bool Regular, IntArray Map );
//! Destructor
~runAnalysis();
//! Run the next analysis
void run(int timestep, std::shared_ptr<Database> db, TwoPhase &Averages, const double *Phi,
void run( int timestep, std::shared_ptr<Database> db, TwoPhase &Averages, const double *Phi,
double *Pressure, double *Velocity, double *fq, double *Den );
void basic( int timestep, std::shared_ptr<Database> db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den );
void WriteVisData(int timestep, std::shared_ptr<Database> vis_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den);
void basic( int timestep, std::shared_ptr<Database> db, SubPhase &Averages, const double *Phi,
double *Pressure, double *Velocity, double *fq, double *Den );
void WriteVisData( int timestep, std::shared_ptr<Database> 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<int,3> d_n; // Number of local cells
std::array<int,3> d_N; // NNumber of local cells with ghosts
std::array<int, 3> d_n; // Number of local cells
std::array<int, 3> 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<std::pair<int, IntArray>> d_last_ids;
std::shared_ptr<std::pair<int, IntArray>> d_last_index;
std::shared_ptr<std::vector<BlobIDType>> d_last_id_map;
std::vector<IO::MeshDataStruct> d_meshData;
std::string d_restartFile;
Utilities::MPI d_comm;
@ -114,8 +119,6 @@ private:
// Friends
friend commWrapper::~commWrapper();
};
#endif

7
common/Communication.h
View File

@ -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<TYPE>& array, int i, int j, int k, TYPE *buffer );
void unpack( Array<TYPE>& array, int i, int j, int k, const TYPE *buffer );
};

26
common/Domain.cpp
View File

@ -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; kp<nprocz; kp++){
for (int jp=0; jp<nprocy; jp++){
for (int ip=0; ip<nprocx; ip++){
@ -609,6 +605,7 @@ void Domain::Decomp( const std::string& Filename )
Comm.send(loc_id,N,rnk,15);
}
// Write the data for this rank data
char LocalRankFilename[40];
sprintf(LocalRankFilename,"ID.%05i",rnk+rank_offset);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(loc_id,1,(nx+2)*(ny+2)*(nz+2),ID);
@ -616,9 +613,8 @@ void Domain::Decomp( const std::string& Filename )
}
}
}
}
else{
delete [] loc_id;
} else {
// Recieve the subdomain from rank = 0
//printf("Ready to recieve data %i at process %i \n", N,rank);
Comm.recv(id.data(),N,0,15);
@ -645,6 +641,7 @@ void Domain::Decomp( const std::string& Filename )
porosity = sum*iVol_global;
if (rank()==0) printf("Media porosity = %f \n",porosity);
//.........................................................
delete [] SegData;
}
void Domain::AggregateLabels( const std::string& filename ){
@ -669,8 +666,7 @@ void Domain::AggregateLabels( const std::string& filename ){
int local_size = (nx-2)*(ny-2)*(nz-2);
long int full_size = long(full_nx)*long(full_ny)*long(full_nz);
signed char *LocalID;
LocalID = new signed char [local_size];
auto LocalID = new signed char [local_size];
//printf("aggregate labels: local size=%i, global size = %i",local_size, full_size);
// assign the ID for the local sub-region
@ -687,8 +683,7 @@ void Domain::AggregateLabels( const std::string& filename ){
// populate the FullID
if (rank() == 0){
signed char *FullID;
FullID = new signed char [full_size];
auto FullID = new signed char [full_size];
// first handle local ID for rank 0
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
@ -727,6 +722,7 @@ void Domain::AggregateLabels( const std::string& filename ){
FILE *OUTFILE = fopen(filename.c_str(),"wb");
fwrite(FullID,1,full_size,OUTFILE);
fclose(OUTFILE);
delete [] FullID;
}
else{
// send LocalID to rank=0
@ -734,8 +730,8 @@ void Domain::AggregateLabels( const std::string& filename ){
int dstrank = 0;
Comm.send(LocalID,local_size,dstrank,tag);
}
delete [] LocalID;
Comm.barrier();
}
/********************************************************
@ -1232,11 +1228,10 @@ void Domain::ReadFromFile(const std::string& Filename,const std::string& Datatyp
// user needs to modify the input file accordingly before LBPM simulator read
// the input file.
//........................................................................................
int rank_offset = 0;
int RANK = rank();
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz;
int64_t global_Nx,global_Ny,global_Nz;
int64_t i,j,k,n;
int64_t i,j,k;
//TODO These offset we may still need them
int64_t xStart,yStart,zStart;
xStart=yStart=zStart=0;
@ -1393,7 +1388,6 @@ void Domain::AggregateLabels( const std::string& filename, DoubleArray &UserData
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
int n = k*nx*ny+j*nx+i;
double local_id_val = UserData(i,j,k);
LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val;
}

109
common/MPI.cpp
View File

@ -401,7 +401,6 @@ MPI_CLASS::MPI_CLASS()
communicator = MPI_CLASS_COMM_NULL;
d_maxTag = mpi_max_tag;
#endif
d_ranks = nullptr;
d_count = nullptr;
d_manage = false;
comm_rank = 0;
@ -435,8 +434,6 @@ void MPI_CLASS::reset()
++N_MPI_Comm_destroyed;
#endif
}
if ( d_ranks != nullptr )
delete[] d_ranks;
delete d_count;
}
if ( d_currentTag == nullptr ) {
@ -448,7 +445,6 @@ void MPI_CLASS::reset()
}
d_manage = false;
d_count = nullptr;
d_ranks = nullptr;
comm_rank = 0;
comm_size = 1;
d_maxTag = 0;
@ -467,7 +463,6 @@ MPI_CLASS::MPI_CLASS( const MPI_CLASS &comm )
d_manage( comm.d_manage ),
comm_rank( comm.comm_rank ),
comm_size( comm.comm_size ),
d_ranks( comm.d_ranks ),
d_maxTag( comm.d_maxTag ),
d_currentTag( comm.d_currentTag )
{
@ -490,7 +485,6 @@ MPI_CLASS::MPI_CLASS( MPI_CLASS &&rhs ) : MPI_CLASS()
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 );
@ -511,7 +505,6 @@ MPI_CLASS &MPI_CLASS::operator=( const MPI_CLASS &comm )
this->communicator = 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<int> 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<int>( 1, myGlobalRank );
if ( d_ranks == nullptr || communicator == MPI_COMM_NULL )
if ( d_isNull )
return std::vector<int>();
// 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<int>( d_ranks, d_ranks + comm_size );
// Check if we are dealing with a serial or global communicator
if ( comm_size == 1 )
return std::vector<int>( 1, globalRank );
if ( comm_size == globalSize ) {
std::vector<int> 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<int>( 1, 1 );
#endif
}
@ -2806,49 +2791,44 @@ MPI_Request MPI_CLASS::IrecvBytes(
}
/************************************************************************
* sendrecv *
************************************************************************/
#if defined( USE_MPI )
template<>
void MPI_CLASS::sendrecv<char>( const char* sendbuf, int sendcount, int dest, int sendtag,
char* recvbuf, int recvcount, int source, int recvtag ) const
void MPI_CLASS::sendrecv<char>( const char *sendbuf, int sendcount, int dest, int sendtag,
char *recvbuf, int recvcount, int source, int recvtag ) const
{
PROFILE_START( "sendrecv<char>", 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<char>", profile_level );
}
template<>
void MPI_CLASS::sendrecv<int>( const int* sendbuf, int sendcount, int dest, int sendtag,
int* recvbuf, int recvcount, int source, int recvtag ) const
void MPI_CLASS::sendrecv<int>( const int *sendbuf, int sendcount, int dest, int sendtag,
int *recvbuf, int recvcount, int source, int recvtag ) const
{
PROFILE_START( "sendrecv<int>", 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<int>", profile_level );
}
template<>
void MPI_CLASS::sendrecv<float>( const float* sendbuf, int sendcount, int dest, int sendtag,
float* recvbuf, int recvcount, int source, int recvtag ) const
void MPI_CLASS::sendrecv<float>( const float *sendbuf, int sendcount, int dest, int sendtag,
float *recvbuf, int recvcount, int source, int recvtag ) const
{
PROFILE_START( "sendrecv<float>", 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<float>", profile_level );
}
template<>
void MPI_CLASS::sendrecv<double>( const double* sendbuf, int sendcount, int dest, int sendtag,
double* recvbuf, int recvcount, int source, int recvtag ) const
void MPI_CLASS::sendrecv<double>( const double *sendbuf, int sendcount, int dest, int sendtag,
double *recvbuf, int recvcount, int source, int recvtag ) const
{
PROFILE_START( "sendrecv<double>", 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<double>", profile_level );
}
#endif
@ -3815,17 +3795,16 @@ MPI MPI::loadBalance( double local, std::vector<double> work )
MPI_ASSERT( (int) work.size() == getSize() );
auto perf = allGather( local );
std::vector<int> I( work.size() );
for ( size_t i=0; i<work.size(); i++)
for ( size_t i = 0; i < work.size(); i++ )
I[i] = i;
auto J = I;
quicksort( perf, I );
quicksort( work, J );
std::vector<int> key( work.size() );
for ( size_t i=0; i<work.size(); i++)
for ( size_t i = 0; i < work.size(); i++ )
key[J[i]] = I[i];
return split( 0, key[getRank()] );
}
} // namespace Utilities

20
common/MPI.h
View File

@ -8,10 +8,14 @@ Copyright (c) 2012 UT-Battelle, LLC
All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
Collection of administrative costs for redistribution of the source code or binary form is allowed. However, collection of a royalty or other fee in excess of good faith amount for cost recovery for such redistribution is prohibited.
Redistribution and use in source and binary forms, with or without modification, are permitted
provided that the following conditions are met: Redistributions of source code must retain the above
copyright notice, this list of conditions and the following disclaimer. Redistributions in binary
form must reproduce the above copyright notice, this list of conditions and the following disclaimer
in the documentation and/or other materials provided with the distribution. Collection of
administrative costs for redistribution of the source code or binary form is allowed. However,
collection of a royalty or other fee in excess of good faith amount for cost recovery for such
redistribution is prohibited.
*/
@ -274,6 +278,8 @@ public: // Member functions
* \brief Return the global ranks for the comm
* \details This returns a vector which contains the global ranks for each
* member of the communicator. The global ranks are defined according to WORLD comm.
* Note: this function is a blocking collective on the current communicator
* (unless the current communicator is global, self, or null)
*/
std::vector<int> globalRanks() const;
@ -796,7 +802,8 @@ public: // Member functions
* @brief This function sends and recieves data using a blocking call
*/
template<class type>
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;

98
common/ScaLBL.cpp
View File

@ -306,9 +306,99 @@ ScaLBL_Communicator::ScaLBL_Communicator(std::shared_ptr <Domain> 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*/
@ -394,7 +484,7 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis
int idx,i,j,k,n;
// Check that Map has size matching sub-domain
if (Map.size(0) != Nx)
if ( (int) Map.size(0) != Nx)
ERROR("ScaLBL_Communicator::MemoryOptimizedLayout: Map array dimensions do not match! \n");
// Initialize Map

31
models/ColorModel.cpp
View File

@ -10,14 +10,31 @@ color lattice boltzmann model
#include <time.h>
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)
@ -408,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;
}
/********************************************************
@ -1097,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();
@ -1334,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();

9
tests/TestNetcdf.cpp
View File

@ -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<int> 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<float> tmp = netcdf::getVar<float>( fid, "tmp" );

102
tests/lbpm_color_simulator.cpp
View File

@ -1,13 +1,13 @@
#include <exception>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#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<Database>( argv[1] );
// Initialize MPI and error handlers
auto multiple = db->getWithDefault<bool>( "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<Database>( 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<std::string>( "TimerFile", "lbpm_color_simulator" );
auto level = db->getWithDefault<int>( "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<std::string>( "TimerFile", "lbpm_color_simulator" );
auto level = db->getWithDefault<int>( "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;
}

1
tests/lbpm_morph_pp.cpp
View File

@ -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();

28
tests/lbpm_uCT_pp.cpp
View File

@ -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<N_levels; i++) {
Array<float> 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<Minkowski>(Dm[0]);
Array <char> phase_label(Nx[0]+2,Ny[0]+2,Nz[0]+2);