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

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

Browse Source
This commit is contained in:
James E McClure
2016-05-27 15:57:33 -04:00
parent 5dcf4f4b4d a77fe2b914
commit 351e7fa593
35 changed files with 2328 additions and 998 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
IO/Reader.cpp
View File

@@ -176,7 +176,7 @@ std::shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const st
var->name = variable;
var->data.resize(N);
if ( precision=="double" ) {
memcpy(var->data.get(),data,bytes);
memcpy(var->data.data(),data,bytes);
} else {
ERROR("Format not implimented");
}

2
IO/Writer.cpp
View File

@@ -120,7 +120,7 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
int dim = mesh.vars[i]->dim;
int type = static_cast<int>(mesh.vars[i]->type);
size_t N = mesh.vars[i]->data.length();
const double* data = N==0 ? NULL:mesh.vars[i]->data.get();
const double* data = N==0 ? NULL:mesh.vars[i]->data.data();
if ( type == static_cast<int>(IO::NullVariable) ) {
ERROR("Variable type not set");
}

93
IO/netcdf.cpp
View File

@@ -32,6 +32,44 @@ inline std::vector<TYPE> reverse( const std::vector<TYPE>& x )
y[i] = x[x.size()-i-1];
return y;
}
// Function to reverse an array
template<class TYPE1, class TYPE2>
inline std::vector<TYPE2> convert( const std::vector<TYPE1>& x )
{
std::vector<TYPE2> y(x.size());
for (size_t i=0; i<x.size(); i++)
y[i] = static_cast<TYPE2>(x[i]);
return y;
}
/****************************************************
* Convert the VariableType to a string *
****************************************************/
std::string VariableTypeName( VariableType type )
{
if ( type == BYTE )
return "BYTE";
else if ( type == SHORT )
return "SHORT";
else if ( type == USHORT )
return "USHORT";
else if ( type == INT )
return "INT";
else if ( type == UINT )
return "UINT";
else if ( type == INT64 )
return "INT64";
else if ( type == UINT64 )
return "UINT64";
else if ( type == FLOAT )
return "FLOAT";
else if ( type == DOUBLE )
return "DOUBLE";
else if ( type == STRING )
return "STRING";
return "Unknown";
}
/****************************************************
@@ -170,70 +208,70 @@ Array<unsigned short> getVar<unsigned short>( int fid, const std::string& var )
{
PROFILE_START("getVar<unsigned short>");
Array<unsigned short> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_ushort( fid, getVarID(fid,var), x.get() );
int err = nc_get_var_ushort( fid, getVarID(fid,var), x.data() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<unsigned short>");
return x;
return x.reverseDim();
}
template<>
Array<short> getVar<short>( int fid, const std::string& var )
{
PROFILE_START("getVar<short>");
Array<short> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_short( fid, getVarID(fid,var), x.get() );
int err = nc_get_var_short( fid, getVarID(fid,var), x.data() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<short>");
return x;
return x.reverseDim();
}
template<>
Array<unsigned int> getVar<unsigned int>( int fid, const std::string& var )
{
PROFILE_START("getVar<unsigned int>");
Array<unsigned int> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_uint( fid, getVarID(fid,var), x.get() );
int err = nc_get_var_uint( fid, getVarID(fid,var), x.data() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<unsigned int>");
return x;
return x.reverseDim();
}
template<>
Array<int> getVar<int>( int fid, const std::string& var )
{
PROFILE_START("getVar<int>");
Array<int> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_int( fid, getVarID(fid,var), x.get() );
int err = nc_get_var_int( fid, getVarID(fid,var), x.data() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<int>");
return x;
return x.reverseDim();
}
template<>
Array<float> getVar<float>( int fid, const std::string& var )
{
PROFILE_START("getVar<float>");
Array<float> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_float( fid, getVarID(fid,var), x.get() );
int err = nc_get_var_float( fid, getVarID(fid,var), x.data() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<float>");
return x;
return x.reverseDim();
}
template<>
Array<double> getVar<double>( int fid, const std::string& var )
{
PROFILE_START("getVar<double>");
Array<double> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_double( fid, getVarID(fid,var), x.get() );
int err = nc_get_var_double( fid, getVarID(fid,var), x.data() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<double>");
return x;
return x.reverseDim();
}
template<>
Array<char> getVar<char>( int fid, const std::string& var )
{
PROFILE_START("getVar<char>");
Array<char> x( reverse(getVarDim(fid,var)) );
int err = nc_get_var_text( fid, getVarID(fid,var), x.get() );
int err = nc_get_var_text( fid, getVarID(fid,var), x.data() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<char>");
return x;
return x.reverseDim();
}
template<>
Array<std::string> getVar<std::string>( int fid, const std::string& var )
@@ -251,6 +289,31 @@ Array<std::string> getVar<std::string>( int fid, const std::string& var )
PROFILE_STOP("getVar<std::string>");
return text;
}
static inline void get_stride_args( const std::vector<int>& start,
const std::vector<int>& count, const std::vector<int>& stride,
size_t *startp, size_t *countp, ptrdiff_t *stridep )
{
for (size_t i=0; i<start.size(); i++)
startp[i] = start[i];
for (size_t i=0; i<count.size(); i++)
countp[i] = count[i];
for (size_t i=0; i<stride.size(); i++)
stridep[i] = stride[i];
}
template<>
Array<short> getVar<short>( int fid, const std::string& var, const std::vector<int>& start,
const std::vector<int>& count, const std::vector<int>& stride )
{
PROFILE_START("getVar<short> (strided)");
Array<short> x( reverse(convert<int,size_t>(count)) );
size_t startp[10], countp[10];
ptrdiff_t stridep[10];
get_stride_args( start, count, stride, startp, countp, stridep );
int err = nc_get_vars_short( fid, getVarID(fid,var), startp, countp, stridep, x.data() );
CHECK_NC_ERR( err );
PROFILE_STOP("getVar<short> (strided)");
return x.reverseDim();
}
/****************************************************
@@ -261,7 +324,7 @@ Array<double> getAtt<double>( int fid, const std::string& att )
{
PROFILE_START("getAtt<double>");
Array<double> x( getAttDim(fid,att) );
int err = nc_get_att_double( fid, NC_GLOBAL, att.c_str(), x.get() );
int err = nc_get_att_double( fid, NC_GLOBAL, att.c_str(), x.data() );
CHECK_NC_ERR( err );
PROFILE_STOP("getAtt<double>");
return x;

17
IO/netcdf.h
View File

@@ -13,6 +13,9 @@ namespace netcdf {
//! Enum to hold variable type
enum VariableType { BYTE, SHORT, USHORT, INT, UINT, INT64, UINT64, FLOAT, DOUBLE, STRING };
//! Convert the VariableType to a string
std::string VariableTypeName( VariableType type );
/*!
* @brief Open netcdf file
@@ -84,6 +87,20 @@ template<class TYPE>
Array<TYPE> getVar( int fid, const std::string& var );
/*!
* @brief Read a strided variable
* @detailed This function reads a strided variable with the given name from the file
* @param fid Handle to the open file
* @param var Variable to read
* @param start Starting corner for the read
* @param count Number of elements to read
* @param stride Stride size for the read
*/
template<class TYPE>
Array<TYPE> getVar( int fid, const std::string& var, const std::vector<int>& start,
const std::vector<int>& count, const std::vector<int>& stride );
/*!
* @brief Read an attribute
* @detailed This function reads an attribute with the given name from the file

8
analysis/analysis.cpp
View File

@@ -52,8 +52,8 @@ int ComputeBlob( const Array<bool>& isPhase, BlobIDArray& LocalBlobID, bool peri
int last = start_id-1;
std::vector<int> neighbor_ids;
neighbor_ids.reserve(N_neighbors);
const bool *isPhasePtr = isPhase.get();
BlobIDType *LocalBlobIDPtr = LocalBlobID.get();
const bool *isPhasePtr = isPhase.data();
BlobIDType *LocalBlobIDPtr = LocalBlobID.data();
for (int z=0; z<Nz; z++) {
for (int y=0; y<Ny; y++) {
for (int x=0; x<Nx; x++) {
@@ -143,7 +143,7 @@ int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
// Compute the local blob ids
size_t N = Nx*Ny*Nz;
Array<bool> isPhase(Nx,Ny,Nz);
memset(isPhase.get(),0,Nx*Ny*Nz*sizeof(bool));
memset(isPhase.data(),0,Nx*Ny*Nz*sizeof(bool));
for (size_t i=0; i<N; i++) {
if ( SignDist(i) <= vS) {
// Solid phase
@@ -765,7 +765,7 @@ void getNewIDs( ID_map_struct& map, BlobIDType& id_max, std::vector<BlobIDType>&
void renumberIDs( const std::vector<BlobIDType>& new_ids, BlobIDArray& IDs )
{
size_t N = IDs.length();
BlobIDType* ids = IDs.get();
BlobIDType* ids = IDs.data();
for (size_t i=0; i<N; i++) {
BlobIDType id = ids[i];
if ( id>=0 )

20
common/Array.cpp
View File

@@ -1,20 +0,0 @@
#include "Array.h"
#include <stdint.h>
/********************************************************
* std::swap *
********************************************************/
namespace std
{
template<> void swap( Array<bool>& v1, Array<bool>& v2 ) { v1.swap(v2); }
template<> void swap( Array<char>& v1, Array<char>& v2 ) { v1.swap(v2); }
template<> void swap( Array<int>& v1, Array<int>& v2 ) { v1.swap(v2); }
template<> void swap( Array<unsigned int>& v1, Array<unsigned int>& v2 ) { v1.swap(v2); }
template<> void swap( Array<int64_t>& v1, Array<int64_t>& v2 ) { v1.swap(v2); }
template<> void swap( Array<uint64_t>& v1, Array<uint64_t>& v2 ) { v1.swap(v2); }
template<> void swap( Array<float>& v1, Array<float>& v2 ) { v1.swap(v2); }
template<> void swap( Array<double>& v1, Array<double>& v2 ) { v1.swap(v2); }
}

404
common/Array.h
View File

@@ -1,48 +1,68 @@
#ifndef included_ArrayClass
#define included_ArrayClass
#include <iostream>
#include <memory>
#include <vector>
#include "shared_ptr.h"
#include "common/Utilities.h"
#include <array>
#include <functional>
#include <iostream>
#include <stdexcept>
#include <memory>
#include <iostream>
#define GET_ARRAY_INDEX(i1,i2,i3,i4) i1+d_N[0]*(i2+d_N[1]*(i3+d_N[2]*i4))
#if defined(DEBUG) || defined(_DEBUG)
#define CHECK_ARRAY_INDEX(i1,i2,i3,i4) \
if ( GET_ARRAY_INDEX(i1,i2,i3,i4)>d_length ) \
ERROR("Index exceeds array bounds");
#define ARRAY_NDIM_MAX 5 // Maximum number of dimensions supported
#define GET_ARRAY_INDEX3D( N, i1, i2, i3 ) i1 + N[0] * ( i2 + N[1] * i3 )
#define GET_ARRAY_INDEX4D( N, i1, i2, i3, i4 ) i1 + N[0] * ( i2 + N[1] * ( i3 + N[2] * i4 ) )
#define GET_ARRAY_INDEX5D( N, i1, i2, i3, i4, i5 ) i1 + N[0] * ( i2 + N[1] * ( i3 + N[2] * ( i4 + N[3] * i5 ) ) )
#if defined( DEBUG ) || defined( _DEBUG )
#define CHECK_ARRAY_INDEX3D( N, i1, i2, i3 ) \
if ( GET_ARRAY_INDEX3D( N, i1, i2, i3 ) < 0 || GET_ARRAY_INDEX3D( N, i1, i2, i3 ) >= d_length ) \
throw std::logic_error( "Index exceeds array bounds" );
#define CHECK_ARRAY_INDEX4D( N, i1, i2, i3, i4 ) \
if ( GET_ARRAY_INDEX4D( N, i1, i2, i3, i4 ) < 0 || \
GET_ARRAY_INDEX4D( N, i1, i2, i3, i4 ) >= d_length ) \
throw std::logic_error( "Index exceeds array bounds" );
#else
#define CHECK_ARRAY_INDEX(i1,i2,i3,i4)
#define CHECK_ARRAY_INDEX3D( N, i1, i2, i3 )
#define CHECK_ARRAY_INDEX4D( N, i1, i2, i3, i4 )
#endif
#if defined( __CUDA_ARCH__ )
#include <cuda.h>
#define HOST_DEVICE __host__ __device__
#else
#define HOST_DEVICE
#endif
/*!
* Class Array is a simple array class
* Class Array is a multi-dimensional array class written by Mark Berrill
*/
template<class TYPE>
template <class TYPE>
class Array
{
public:
/*!
* Create a new empty Array
*/
Array( );
Array();
/*!
* Create a new 1D Array with the given number of elements
* @param N Number of elements in the array
*/
Array( size_t N );
explicit Array( size_t N );
/*!
* Create a new 2D Array with the given number of rows and columns
* @param N_rows Number of rows
* @param N_columns Number of columns
*/
Array( size_t N_rows, size_t N_columns );
explicit Array( size_t N_rows, size_t N_columns );
/*!
* Create a new 3D Array with the given number of rows and columns
@@ -50,33 +70,52 @@ public:
* @param N2 Number of columns
* @param N3 Number of elements in the third dimension
*/
Array( size_t N1, size_t N2, size_t N3 );
explicit Array( size_t N1, size_t N2, size_t N3 );
/*!
* Create a multi-dimensional Array with the given number of elements
* @param N Number of elements in each dimension
* @param data Optional raw array to copy the src data
*/
Array( const std::vector<size_t>& N );
explicit Array( const std::vector<size_t> &N, const TYPE *data = NULL );
/*!
* Copy constructor
* @param rhs Array to copy
*/
Array( const Array& rhs );
Array( const Array &rhs );
/*!
* Move constructor
* @param rhs Array to copy
*/
Array( Array &&rhs );
/*!
* Assignment operator
* @param rhs Array to copy
*/
Array& operator=( const Array& rhs );
Array &operator=( const Array &rhs );
/*!
* Move assignment operator
* @param rhs Array to copy
*/
Array &operator=( Array &&rhs );
/*!
* Assignment operator
* @param rhs std::vector to copy
*/
Array &operator=( const std::vector<TYPE> &rhs );
/*!
* Create a 1D Array view to a raw block of data
* @param N Number of elements in the array
* @param data Pointer to the data
*/
static std::shared_ptr<Array> view( size_t N, std::shared_ptr<TYPE> data );
static std::shared_ptr<Array> view( size_t N, std::shared_ptr<TYPE> const &data );
/*!
* Create a new 2D Array with the given number of rows and columns
@@ -84,7 +123,8 @@ public:
* @param N_columns Number of columns
* @param data Pointer to the data
*/
static std::shared_ptr<Array> view( size_t N_rows, size_t N_columns, std::shared_ptr<TYPE> data );
static std::shared_ptr<Array> view(
size_t N_rows, size_t N_columns, std::shared_ptr<TYPE> const &data );
/*!
* Create a new 3D Array view to a raw block of data
@@ -93,22 +133,25 @@ public:
* @param N3 Number of elements in the third dimension
* @param data Pointer to the data
*/
static std::shared_ptr<Array> view( size_t N1, size_t N2, size_t N3, std::shared_ptr<TYPE> data );
static std::shared_ptr<Array> view(
size_t N1, size_t N2, size_t N3, std::shared_ptr<TYPE> const &data );
/*!
* Create a multi-dimensional Array view to a raw block of data
* @param N Number of elements in each dimension
* @param data Pointer to the data
*/
static std::shared_ptr<Array> view( const std::vector<size_t>& N, std::shared_ptr<TYPE> data );
static std::shared_ptr<Array> view(
const std::vector<size_t> &N, std::shared_ptr<TYPE> const &data );
/*!
* Create a 1D Array view to a raw block of data
* @param N Number of elements in the array
* @param data Pointer to the data
*/
static std::shared_ptr<const Array> constView( size_t N, std::shared_ptr<const TYPE> data );
static std::shared_ptr<const Array> constView(
size_t N, std::shared_ptr<const TYPE> const &data );
/*!
* Create a new 2D Array with the given number of rows and columns
@@ -116,7 +159,8 @@ public:
* @param N_columns Number of columns
* @param data Pointer to the data
*/
static std::shared_ptr<const Array> constView( size_t N_rows, size_t N_columns, std::shared_ptr<const TYPE> data );
static std::shared_ptr<const Array> constView(
size_t N_rows, size_t N_columns, std::shared_ptr<const TYPE> const &data );
/*!
* Create a new 3D Array view to a raw block of data
@@ -125,63 +169,121 @@ public:
* @param N3 Number of elements in the third dimension
* @param data Pointer to the data
*/
static std::shared_ptr<const Array> constView( size_t N1, size_t N2, size_t N3, std::shared_ptr<const TYPE> data );
static std::shared_ptr<const Array> constView(
size_t N1, size_t N2, size_t N3, std::shared_ptr<const TYPE> const &data );
/*!
* Create a multi-dimensional Array view to a raw block of data
* @param N Number of elements in each dimension
* @param data Pointer to the data
*/
static std::shared_ptr<const Array> constView( const std::vector<size_t>& N, std::shared_ptr<const TYPE> data );
static std::shared_ptr<const Array> constView(
const std::vector<size_t> &N, std::shared_ptr<const TYPE> const &data );
/*!
* Make this object a view of the src
* @param src Source vector to take the view of
*/
void view2( Array &src );
/*!
* Make this object a view of the data
* @param N Number of elements in each dimension
* @param data Pointer to the data
*/
void view2( const std::vector<size_t> &N, std::shared_ptr<TYPE> const &data );
/*!
* Make this object a view of the raw data (expert use only).
* Use view2( N, std::shared_ptr(data,[](TYPE*){}) ) instead.
* Note: this interface is not recommended as it does not protect from
* the src data being deleted while still being used by the Array.
* Additionally for maximum performance it does not set the internal shared_ptr
* so functions like getPtr and resize will not work correctly.
* @param N Number of elements in each dimension
* @param data Pointer to the data
*/
void viewRaw( const std::initializer_list<size_t> &N, TYPE *data );
/*!
* Make this object a view of the raw data (expert use only).
* Use view2( N, std::shared_ptr(data,[](TYPE*){}) ) instead.
* Note: this interface is not recommended as it does not protect from
* the src data being deleted while still being used by the Array.
* Additionally for maximum performance it does not set the internal shared_ptr
* so functions like getPtr and resize will not work correctly.
* @param N Number of elements in each dimension
* @param data Pointer to the data
*/
void viewRaw( const std::vector<size_t> &N, TYPE *data );
/*!
* Convert an array of one type to another. This may or may not allocate new memory.
* @param array Input array
*/
template <class TYPE2>
static std::shared_ptr<Array<TYPE2>> convert( std::shared_ptr<Array<TYPE>> array );
/*!
* Convert an array of one type to another. This may or may not allocate new memory.
* @param array Input array
*/
template<class TYPE2>
static std::shared_ptr<Array<TYPE2> > convert( std::shared_ptr<Array<TYPE> > array );
/*!
* Convert an array of one type to another. This may or may not allocate new memory.
* @param array Input array
*/
template<class TYPE2>
static std::shared_ptr<const Array<TYPE2> > convert( std::shared_ptr<const Array<TYPE> > array );
template <class TYPE2>
static std::shared_ptr<const Array<TYPE2>> convert( std::shared_ptr<const Array<TYPE>> array );
/*!
* Copy and convert data from another array to this array
* @param array Source array
*/
template<class TYPE2>
void copy( const Array<TYPE2>& array );
template <class TYPE2>
void copy( const Array<TYPE2> &array );
/*!
* Copy and convert data from a raw vector to this array.
* Note: The current array must be allocated to the proper size first.
* @param array Source array
*/
template<class TYPE2>
void copy( const TYPE2* array );
template <class TYPE2>
void copy( const TYPE2 *array );
/*!
* Copy and convert data from this array to a raw vector.
* @param array Source array
*/
template <class TYPE2>
void copyTo( TYPE2 *array ) const;
/*!
* Fill the array with the given value
* @param value Value to fill
*/
void fill( const TYPE& value );
void fill( const TYPE &value );
/*!
* Scale the array by the given value
* @param scale Value to scale by
*/
void scale( const TYPE &scale );
//! Destructor
~Array( );
~Array();
//! Clear the data in the array
void clear();
//! Return the size of the Array
inline int ndim( ) const { return d_ndim; }
inline int ndim() const { return d_ndim; }
//! Return the size of the Array
inline std::vector<size_t> size( ) const { return std::vector<size_t>(d_N,d_N+d_ndim); }
std::vector<size_t> size() const;
//! Return the size of the Array
@@ -189,11 +291,11 @@ public:
//! Return the size of the Array
inline size_t length( ) const { return d_length; }
inline size_t length() const { return d_length; }
//! Return true if the Array is empty
inline bool empty( ) const { return d_length==0; }
inline bool empty() const { return d_length == 0; }
/*!
@@ -221,126 +323,232 @@ public:
* Resize the Array
* @param N Number of elements in each dimension
*/
void resize( const std::vector<size_t>& N );
void resize( const std::vector<size_t> &N );
/*!
* Resize the given dimension of the array
* @param dim The dimension to resize
* @param N Number of elements for the given dimension
* @param value Value to initialize new elements
*/
void resizeDim( int dim, size_t N, const TYPE &value );
/*!
* Reshape the Array (total size of array will not change)
* @param N Number of elements in each dimension
*/
void reshape( const std::vector<size_t>& N );
void reshape( const std::vector<size_t> &N );
/*!
* Subset the Array (total size of array will not change)
* @param index Index to subset (imin,imax,jmin,jmax,kmin,kmax,...)
*/
template<class TYPE2=TYPE>
Array<TYPE2> subset( const std::vector<size_t> &index ) const;
/*!
* Copy data from an array into a subset of this array
* @param index Index of the subset (imin,imax,jmin,jmax,kmin,kmax,...)
* @param subset The subset array to copy from
*/
template <class TYPE2>
void copySubset( const std::vector<size_t> &index, const Array<TYPE2> &subset );
/*!
* Access the desired element
* @param i The row index
*/
HOST_DEVICE inline TYPE &operator()( size_t i )
{
CHECK_ARRAY_INDEX3D( d_N, i, 0, 0 ) return d_data[i];
}
/*!
* Access the desired element
* @param i The row index
*/
HOST_DEVICE inline const TYPE &operator()( size_t i ) const
{
CHECK_ARRAY_INDEX3D( d_N, i, 0, 0 ) return d_data[i];
}
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline TYPE& operator()( size_t i ) { CHECK_ARRAY_INDEX(i,0,0,0) return d_data[i]; }
HOST_DEVICE inline TYPE &operator()( size_t i, size_t j )
{
CHECK_ARRAY_INDEX3D( d_N, i, j, 0 ) return d_data[i + j * d_N[0]];
}
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline const TYPE& operator()( size_t i ) const { CHECK_ARRAY_INDEX(i,0,0,0) return d_data[i]; }
HOST_DEVICE inline const TYPE &operator()( size_t i, size_t j ) const
{
CHECK_ARRAY_INDEX3D( d_N, i, j, 0 ) return d_data[i + j * d_N[0]];
}
/*!
* Access the desired element
* @param i The row index
* @param j The column index
* @param k The third index
*/
inline TYPE& operator()( size_t i, size_t j ) { CHECK_ARRAY_INDEX(i,j,0,0) return d_data[i+j*d_N[0]]; }
HOST_DEVICE inline TYPE &operator()( size_t i, size_t j, size_t k )
{
CHECK_ARRAY_INDEX3D( d_N, i, j, k ) return d_data[GET_ARRAY_INDEX3D( d_N, i, j, k )];
}
/*!
* Access the desired element
* @param i The row index
* @param j The column index
* @param k The third index
*/
inline const TYPE& operator()( size_t i, size_t j ) const { CHECK_ARRAY_INDEX(i,j,0,0) return d_data[i+j*d_N[0]]; }
HOST_DEVICE inline const TYPE &operator()( size_t i, size_t j, size_t k ) const
{
CHECK_ARRAY_INDEX3D( d_N, i, j, k ) return d_data[GET_ARRAY_INDEX3D( d_N, i, j, k )];
}
/*!
* Access the desired element
* @param i The row index
* @param j The column index
* @param k The third index
* @param l The fourth index
*/
inline TYPE& operator()( size_t i, size_t j, size_t k ) { CHECK_ARRAY_INDEX(i,j,k,0) return d_data[GET_ARRAY_INDEX(i,j,k,0)]; }
HOST_DEVICE inline TYPE &operator()( size_t i, size_t j, size_t k, size_t l )
{
CHECK_ARRAY_INDEX4D( d_N, i, j, k, l ) return d_data[GET_ARRAY_INDEX4D( d_N, i, j, k, l )];
}
/*!
* Access the desired element
* @param i The row index
* @param j The column index
* @param k The third index
* @param l The fourth index
*/
inline const TYPE& operator()( size_t i, size_t j, size_t k ) const { CHECK_ARRAY_INDEX(i,j,k,0) return d_data[GET_ARRAY_INDEX(i,j,k,0)]; }
HOST_DEVICE inline const TYPE &operator()( size_t i, size_t j, size_t k, size_t l ) const
{
CHECK_ARRAY_INDEX4D( d_N, i, j, k, l ) return d_data[GET_ARRAY_INDEX4D( d_N, i, j, k, l )];
}
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline TYPE& operator()( size_t i, size_t j, size_t k, size_t m ) { CHECK_ARRAY_INDEX(i,j,k,m) return d_data[GET_ARRAY_INDEX(i,j,k,m)]; }
//! Check if two matrices are equal
// Equality means the dimensions and data have to be identical
bool operator==( const Array &rhs ) const;
/*!
* Access the desired element
* @param i The row index
* @param j The column index
*/
inline const TYPE& operator()( size_t i, size_t j, size_t k, size_t m ) const { CHECK_ARRAY_INDEX(i,j,k,m) return d_data[GET_ARRAY_INDEX(i,j,k,m)]; }
//! Check if two matricies are equal
bool operator==( const Array& rhs ) const;
//! Check if two matricies are not equal
inline bool operator!=( const Array& rhs ) const { return !this->operator==(rhs); }
//! Check if two matrices are not equal
inline bool operator!=( const Array &rhs ) const { return !this->operator==( rhs ); }
//! Return the pointer to the raw data
inline std::shared_ptr<TYPE> getPtr( ) { return d_ptr; }
inline std::shared_ptr<TYPE> getPtr() { return d_ptr; }
//! Return the pointer to the raw data
inline std::shared_ptr<const TYPE> getPtr( ) const { return d_ptr; }
inline std::shared_ptr<const TYPE> getPtr() const { return d_ptr; }
//! Return the pointer to the raw data
inline TYPE* get( ) { return d_data; }
HOST_DEVICE inline TYPE *data() { return d_data; }
//! Return the pointer to the raw data
inline const TYPE* get( ) const { return d_data; }
HOST_DEVICE inline const TYPE *data() const { return d_data; }
//! Return true if NaNs are present
inline bool NaNs( ) const;
inline bool NaNs() const;
//! Return the smallest value
inline TYPE min( ) const;
inline TYPE min() const;
//! Return the largest value
inline TYPE max( ) const;
inline TYPE max() const;
//! Return the sum of all elements
inline TYPE sum( ) const;
inline TYPE sum() const;
//! Return the mean of all elements
inline TYPE mean() const;
//! Return the min of all elements in a given direction
Array<TYPE> min( int dir ) const;
//! Return the max of all elements in a given direction
Array<TYPE> max( int dir ) const;
//! Return the sum of all elements in a given direction
std::shared_ptr<Array<TYPE> > sum( int dir ) const;
Array<TYPE> sum( int dir ) const;
//! Swap the data in this with rhs
inline void swap( Array& rhs );
//! Return the smallest value
inline TYPE min( const std::vector<size_t> &index ) const;
//! Return the largest value
inline TYPE max( const std::vector<size_t> &index ) const;
//! Return the sum of all elements
inline TYPE sum( const std::vector<size_t> &index ) const;
//! Return the mean of all elements
inline TYPE mean( const std::vector<size_t> &index ) const;
//! Find all elements that match the operator
std::vector<size_t> find(
const TYPE &value, std::function<bool( const TYPE &, const TYPE & )> compare ) const;
//! Add another array
Array &operator+=( const Array &rhs );
//! Subtract another array
Array &operator-=( const Array &rhs );
//! Add a scalar
Array &operator+=( const TYPE &rhs );
//! Subtract a scalar
Array &operator-=( const TYPE &rhs );
//! Print an array
void print( std::ostream& os, const std::string& name="A", const std::string& prefix="" ) const;
//! Multiply two arrays
static Array multiply( const Array& a, const Array& b );
//! Transpose an array
Array<TYPE> reverseDim( ) const;
//! Coarsen an array using the given filter
Array<TYPE> coarsen( const Array<TYPE>& filter ) const;
private:
int d_ndim;
size_t d_N[4];
size_t d_length;
TYPE *d_data;
std::shared_ptr<TYPE> d_ptr;
void allocate( const std::vector<size_t>& N );
int d_ndim; // Number of dimensions in array
size_t d_N[ARRAY_NDIM_MAX]; // Size of each dimension
size_t d_length; // Total length of array
TYPE *d_data; // Raw pointer to data in array
std::shared_ptr<TYPE> d_ptr; // Shared pointer to data in array
void allocate( const std::vector<size_t> &N );
private:
template<class TYPE2>
inline bool sizeMatch( const Array<TYPE2>& rhs ) const;
inline void checkSubsetIndex( const std::vector<size_t> &index ) const;
inline std::array<size_t, 5> getDimArray() const;
static inline void getSubsetArrays( const std::vector<size_t> &index,
std::array<size_t, 5> &first, std::array<size_t, 5> &last, std::array<size_t, 5> &N );
};
typedef Array<int> IntArray;
typedef Array<double> DoubleArray;
typedef Array<float> FloatArray;
typedef Array<int> IntArray;
#include "common/Array.hpp"
#endif

1026
common/Array.hpp
View File

File diff suppressed because it is too large Load Diff

5
common/Communication.h
View File

@@ -2,7 +2,8 @@
#define COMMUNICATION_H_INC
#include "common/MPI_Helpers.h"
#include "Array.h"
#include "common/Utilities.h"
#include "common/Array.h"
// ********** COMMUNICTION **************************************
/*
@@ -287,7 +288,7 @@ inline void CommunicateMeshHalo(DoubleArray &Mesh, MPI_Comm Communicator,
{
int sendtag, recvtag;
sendtag = recvtag = 7;
double *MeshData = Mesh.get();
double *MeshData = Mesh.data();
PackMeshData(sendList_x, sendCount_x ,sendbuf_x, MeshData);
PackMeshData(sendList_X, sendCount_X ,sendbuf_X, MeshData);
PackMeshData(sendList_y, sendCount_y ,sendbuf_y, MeshData);

1
common/Communication.hpp
View File

@@ -3,6 +3,7 @@
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "common/Utilities.h"
#include "ProfilerApp.h"

4
common/Domain.cpp
View File

@@ -547,7 +547,7 @@ void Domain::CommunicateMeshHalo(DoubleArray &Mesh)
{
int sendtag, recvtag;
sendtag = recvtag = 7;
double *MeshData = Mesh.get();
double *MeshData = Mesh.data();
PackMeshData(sendList_x, sendCount_x ,sendData_x, MeshData);
PackMeshData(sendList_X, sendCount_X ,sendData_X, MeshData);
PackMeshData(sendList_y, sendCount_y ,sendData_y, MeshData);
@@ -631,7 +631,7 @@ void Domain::BlobComm(MPI_Comm Communicator)
int sendtag, recvtag;
sendtag = recvtag = 51;
//......................................................................................
int *BlobLabelData = BlobLabel.get();
int *BlobLabelData = BlobLabel.data();
PackBlobData(sendList_x, sendCount_x ,sendBuf_x, BlobLabelData);
PackBlobData(sendList_X, sendCount_X ,sendBuf_X, BlobLabelData);
PackBlobData(sendList_y, sendCount_y ,sendBuf_y, BlobLabelData);

8
common/TwoPhase.cpp
View File

@@ -790,8 +790,8 @@ void TwoPhase::ComponentAverages()
}
*/
MPI_Barrier(Dm.Comm);
MPI_Allreduce(ComponentAverages_NWP.get(),RecvBuffer.get(),BLOB_AVG_COUNT*NumberComponents_NWP, MPI_DOUBLE,MPI_SUM,Dm.Comm);
// MPI_Reduce(ComponentAverages_NWP.get(),RecvBuffer.get(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
MPI_Allreduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT*NumberComponents_NWP, MPI_DOUBLE,MPI_SUM,Dm.Comm);
// MPI_Reduce(ComponentAverages_NWP.data(),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
if (Dm.rank==0){
printf("rescaling... \n");
@@ -888,8 +888,8 @@ void TwoPhase::ComponentAverages()
// reduce the wetting phase averages
for (int b=0; b<NumberComponents_WP; b++){
MPI_Barrier(Dm.Comm);
// MPI_Allreduce(&ComponentAverages_WP(0,b),RecvBuffer.get(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Reduce(&ComponentAverages_WP(0,b),RecvBuffer.get(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
// MPI_Allreduce(&ComponentAverages_WP(0,b),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Reduce(&ComponentAverages_WP(0,b),RecvBuffer.data(),BLOB_AVG_COUNT,MPI_DOUBLE,MPI_SUM,0,Dm.Comm);
for (int idx=0; idx<BLOB_AVG_COUNT; idx++) ComponentAverages_WP(idx,b)=RecvBuffer(idx);
}

159
common/imfilter.h Normal file
View File

@@ -0,0 +1,159 @@
// These functions mimic the behavior of imfilter in MATLAB
#ifndef included_imfilter
#define included_imfilter
#include "common/Utilities.h"
#include "common/Array.h"
#include <vector>
namespace imfilter {
//! enum to store the BC type
enum class BC { fixed=0, symmetric=1, replicate=2, circular=3 };
/*!
* @brief N-D filtering of multidimensional images
* @details imfilter filters the multidimensional array A with the
* multidimensional filter H. The result B has the same size and class as A.
* @param[in] A The input array (Nx,Ny,Nz)
* @param[in] H The filter (2*Nhx+1,2*Nhy+1,...)
* @param[in] boundary The boundary conditions to apply (ndim):
* fixed - Input array values outside the bounds of the array are
* implicitly assumed to have the value X
* symmetric - Input array values outside the bounds of the array are
* computed by mirror-reflecting the array across the array border
* replicate - Input array values outside the bounds of the array are
* assumed to equal the nearest array border value
* circular - Input array values outside the bounds of the array are
* computed by implicitly assuming the input array is periodic.
* @param[in] X The value to use for boundary conditions (only used if boundary==fixed)
*/
template<class TYPE>
Array<TYPE> imfilter( const Array<TYPE>& A, const Array<TYPE>& H, const std::vector<imfilter::BC>& boundary, const TYPE X=0 );
/*!
* @brief N-D filtering of multidimensional images
* @details imfilter filters the multidimensional array A with the
* multidimensional filter H. The result B has the same size and class as A.
* @param[in] A The input array (Nx,Ny,Nz)
* @param[in] Nh The size of the filter
* @param[in] H The filter function to use ( y = H(data) )
* Note that the data passed to this function will be of size 2*Nh+1
* @param[in] boundary The boundary conditions to apply (ndim):
* fixed - Input array values outside the bounds of the array are
* implicitly assumed to have the value X
* symmetric - Input array values outside the bounds of the array are
* computed by mirror-reflecting the array across the array border
* replicate - Input array values outside the bounds of the array are
* assumed to equal the nearest array border value
* circular - Input array values outside the bounds of the array are
* computed by implicitly assuming the input array is periodic.
* @param[in] X The value to use for boundary conditions (only used if boundary==fixed)
*/
template<class TYPE>
Array<TYPE> imfilter( const Array<TYPE>& A, const std::vector<int>& Nh,
std::function<TYPE(const Array<TYPE>&)> H,
const std::vector<imfilter::BC>& boundary, const TYPE X=0 );
/*!
* @brief N-D filtering of multidimensional images
* @details imfilter filters the multidimensional array A with the
* multidimensional filter H. The result B has the same size and class as A.
* This version works with separable filters and is more efficient than a single filter.
* @param[in] A The input array (Nx,Ny,Nz)
* @param[in] H The filter [2*Nhx+1,2*Nhy+1,...]
* @param[in] boundary The boundary conditions to apply (ndim):
* fixed - Input array values outside the bounds of the array are
* implicitly assumed to have the value X
* symmetric - Input array values outside the bounds of the array are
* computed by mirror-reflecting the array across the array border
* replicate - Input array values outside the bounds of the array are
* assumed to equal the nearest array border value
* circular - Input array values outside the bounds of the array are
* computed by implicitly assuming the input array is periodic.
* @param[in] X The value to use for boundary conditions (only used if boundary==fixed)
*/
template<class TYPE>
Array<TYPE> imfilter_separable( const Array<TYPE>& A, const std::vector<Array<TYPE>>& H,
const std::vector<imfilter::BC>& boundary, const TYPE X=0 );
/*!
* @brief N-D filtering of multidimensional images
* @details imfilter filters the multidimensional array A with the
* multidimensional filter H. The result B has the same size and class as A.
* This version works with separable filters and is more efficient than a single filter.
* @param[in] A The input array (Nx,Ny,Nz)
* @param[in] H The filter [2*Nhx+1,2*Nhy+1,...]
* @param[in] boundary The boundary conditions to apply (ndim):
* fixed - Input array values outside the bounds of the array are
* implicitly assumed to have the value X
* symmetric - Input array values outside the bounds of the array are
* computed by mirror-reflecting the array across the array border
* replicate - Input array values outside the bounds of the array are
* assumed to equal the nearest array border value
* circular - Input array values outside the bounds of the array are
* computed by implicitly assuming the input array is periodic.
* @param[in] X The value to use for boundary conditions (only used if boundary==fixed)
*/
template<class TYPE>
Array<TYPE> imfilter_separable( const Array<TYPE>& A, const std::vector<int>& Nh,
std::vector<std::function<TYPE(const Array<TYPE>&)>> H,
const std::vector<imfilter::BC>& boundary, const TYPE X=0 );
/*!
* @brief N-D filtering of multidimensional images
* @details imfilter filters the multidimensional array A with the
* multidimensional filter H. The result B has the same size and class as A.
* This version works with separable filters and is more efficient than a single filter.
* @param[in] A The input array (Nx,Ny,Nz)
* @param[in] H The filter [2*Nhx+1,2*Nhy+1,...]
* @param[in] boundary The boundary conditions to apply (ndim):
* fixed - Input array values outside the bounds of the array are
* implicitly assumed to have the value X
* symmetric - Input array values outside the bounds of the array are
* computed by mirror-reflecting the array across the array border
* replicate - Input array values outside the bounds of the array are
* assumed to equal the nearest array border value
* circular - Input array values outside the bounds of the array are
* computed by implicitly assuming the input array is periodic.
* @param[in] X The value to use for boundary conditions (only used if boundary==fixed)
*/
template<class TYPE>
Array<TYPE> imfilter_separable( const Array<TYPE>& A, const std::vector<int>& Nh,
std::vector<std::function<TYPE(int, const TYPE*)>> H,
const std::vector<imfilter::BC>& boundary, const TYPE X=0 );
/**
* @brief Create a filter to use with imfilter
* @details This function creates one of several predefined filters
* to use with imfilter. The filter will always sum to 1.
* Note: this function allocates memory with the new command, the user must call delete.
*
* @param[in] N The stencil size in each direction
* @param[in] type The type of filter to create
* average - Simple averaging filter
* gaussian - Gaussian filter with given standard deviation.
* Optional argument is a double array of size ndim
* giving the standard deviation in each direction.
* A default value of 0.5 is used if not provided.
* \param[in] args An optional argument that some of the filters use
*/
template<class TYPE>
Array<TYPE> create_filter( const std::vector<int>& N, const std::string &type, const void *args = NULL );
}
#include "common/imfilter.hpp"
#endif

378
common/imfilter.hpp Normal file
View File

@@ -0,0 +1,378 @@
#include "common/imfilter.h"
#include "ProfilerApp.h"
#include <math.h>
#include <string.h>
#define IMFILTER_INSIST INSIST
#define IMFILTER_ASSERT ASSERT
#define IMFILTER_ERROR ERROR
// Function to convert an index
static inline int imfilter_index( int index, const int N, const imfilter::BC bc )
{
if ( index < 0 || index >= N ) {
if ( bc == imfilter::BC::symmetric ) {
index = ( 2 * N - index ) % N;
} else if ( bc == imfilter::BC::replicate ) {
index = index < 0 ? 0 : N - 1;
} else if ( bc == imfilter::BC::circular ) {
index = ( index + N ) % N;
} else if ( bc == imfilter::BC::fixed ) {
index = -1;
}
}
return index;
}
// Function to copy a 1D array and pad with the appropriate BC
template<class TYPE>
static inline void copy_array( const int N, const int Ns, const int Nh,
const TYPE *A, const imfilter::BC BC, const TYPE X, TYPE *B )
{
// Fill the center with a memcpy
for (int i=0; i<N; i++ )
B[i+Nh] = A[i*Ns];
// Fill the boundaries
for (int i=0; i<Nh; i++ ) {
int j1 = imfilter_index( -(i+1), N, BC );
int j2 = imfilter_index( N+i, N, BC );
B[Nh-i-1] = j1==-1 ? X : B[Nh+j1];
B[N+Nh+i] = j2==-1 ? X : B[Nh+j2];
}
}
/********************************************************
* Perform a 1D filter in a single direction *
********************************************************/
template<class TYPE>
static void filter_direction( int Ns, int N, int Ne, int Nh, const TYPE *H,
imfilter::BC boundary, TYPE X, TYPE *A )
{
if ( Nh < 0 )
IMFILTER_ERROR("Invalid filter size");
if ( Nh == 0 ) {
for (int i=0; i<Ns*N*Ne; i++)
A[i] *= H[0];
return;
}
TYPE *tmp = new TYPE[N+2*Nh];
for (int j=0; j<Ne; j++) {
for (int i=0; i<Ns; i++) {
copy_array( N, Ns, Nh, &A[i+j*Ns*N], boundary, X, tmp );
for (int k=0; k<N; k++) {
TYPE tmp2 = 0;
for (int m=0; m<=2*Nh; m++)
tmp2 += H[m] * tmp[k+m];
A[i+k*Ns+j*Ns*N] = tmp2;
}
}
}
delete[] tmp;
}
template<class TYPE>
static void filter_direction( int Ns, int N, int Ne, int Nh,
std::function<TYPE(const Array<TYPE>&)> H, imfilter::BC boundary, TYPE X, TYPE *A )
{
if ( Nh < 0 )
IMFILTER_ERROR("Invalid filter size");
TYPE *tmp = new TYPE[N+2*Nh];
Array<TYPE> tmp2(2*Nh+1);
for (int j=0; j<Ne; j++) {
for (int i=0; i<Ns; i++) {
copy_array( N, Ns, Nh, &A[i+j*Ns*N], boundary, X, tmp );
for (int k=0; k<N; k++) {
for (int m=0; m<=2*Nh; m++)
tmp2(m) = tmp[k+m];
A[i+k*Ns+j*Ns*N] = H(tmp2);
}
}
}
delete[] tmp;
}
template<class TYPE>
static void filter_direction( int Ns, int N, int Ne, int Nh,
std::function<TYPE(int, const TYPE*)> H, imfilter::BC boundary, TYPE X, TYPE *A )
{
if ( Nh < 0 )
IMFILTER_ERROR("Invalid filter size");
TYPE *tmp = new TYPE[N+2*Nh];
int Nh2 = 2*Nh+1;
for (int j=0; j<Ne; j++) {
for (int i=0; i<Ns; i++) {
copy_array( N, Ns, Nh, &A[i+j*Ns*N], boundary, X, tmp );
for (int k=0; k<N; k++)
A[i+k*Ns+j*Ns*N] = H(Nh2,&tmp[k]);
}
}
delete[] tmp;
}
/********************************************************
* Create a filter *
********************************************************/
template<class TYPE>
Array<TYPE> imfilter::create_filter( const std::vector<int>& N0, const std::string &type, const void *args )
{
std::vector<size_t> N2(N0.size());
for (size_t i=0; i<N2.size(); i++)
N2[i] = 2*N0[i]+1;
Array<TYPE> h(N2);
h.fill(0);
if ( type == "average" ) {
// average
h.fill( 1.0 / static_cast<TYPE>( h.length() ) );
} else if ( type == "gaussian" ) {
// gaussian
if ( N0.size() > 3 )
IMFILTER_ERROR( "Not implimented for dimensions > 3" );
TYPE std[3] = { 0.5, 0.5, 0.5 };
if ( args != NULL ) {
const TYPE *args2 = reinterpret_cast<const TYPE*>( args );
for ( int d = 0; d < N0.size(); d++ )
std[d] = args2[d];
}
auto N = N0;
N.resize(3,0);
for ( int k = -N[2]; k <= N[2]; k++ ) {
for ( int j = -N[1]; j <= N[1]; j++ ) {
for ( int i = -N[0]; i <= N[0]; i++ ) {
h(i+N[0],j+N[1],k+N[2]) =
exp( -i * i / ( 2 * std[0] * std[0] ) ) *
exp( -j * j / ( 2 * std[1] * std[1] ) ) *
exp( -k * k / ( 2 * std[2] * std[2] ) );
}
}
}
h.scale( 1.0/h.sum() );
} else {
IMFILTER_ERROR( "Unknown filter" );
}
return h;
}
// Perform 2-D filtering
template<class TYPE>
void imfilter_2D( int Nx, int Ny, const TYPE *A, int Nhx, int Nhy, const TYPE *H,
imfilter::BC BCx, imfilter::BC BCy, const TYPE X, TYPE *B )
{
IMFILTER_ASSERT( A != B );
PROFILE_START( "imfilter_2D" );
memset( B, 0, Nx * Ny * sizeof( TYPE ) );
for ( int j1 = 0; j1 < Ny; j1++ ) {
for ( int i1 = 0; i1 < Nx; i1++ ) {
TYPE tmp = 0;
if ( i1 >= Nhx && i1 < Nx - Nhx && j1 >= Nhy && j1 < Ny - Nhy ) {
int ijkh = 0;
for ( int j2 = j1 - Nhy; j2 <= j1 + Nhy; j2++ ) {
for ( int i2 = i1 - Nhx; i2 <= i1 + Nhx; i2++, ijkh++ )
tmp += H[ijkh] * A[i2 + j2 * Nx];
}
} else {
int ijkh = 0;
for ( int jh = -Nhy; jh <= Nhy; jh++ ) {
int j2 = imfilter_index( j1+jh, Ny, BCy );
for ( int ih = -Nhx; ih <= Nhx; ih++ ) {
int i2 = imfilter_index( i1+ih, Nx, BCx );
bool fixed = i2 == -1 || j2 == -1;
TYPE A2 = fixed ? X : A[i2 + j2 * Nx];
tmp += H[ijkh] * A2;
ijkh++;
}
}
}
B[i1 + j1 * Nx] = tmp;
}
}
PROFILE_STOP( "imfilter_2D" );
}
// Perform 3-D filtering
template<class TYPE>
void imfilter_3D( int Nx, int Ny, int Nz, const TYPE *A, int Nhx, int Nhy, int Nhz,
const TYPE *H, imfilter::BC BCx, imfilter::BC BCy, imfilter::BC BCz,
const TYPE X, TYPE *B )
{
IMFILTER_ASSERT( A != B );
PROFILE_START( "imfilter_3D" );
memset( B, 0, Nx * Ny * Nz * sizeof( TYPE ) );
for ( int k1 = 0; k1 < Nz; k1++ ) {
for ( int j1 = 0; j1 < Ny; j1++ ) {
for ( int i1 = 0; i1 < Nx; i1++ ) {
TYPE tmp = 0;
int ijkh = 0;
for ( int kh = -Nhz; kh <= Nhz; kh++ ) {
int k2 = imfilter_index( k1+kh, Nz, BCz );
for ( int jh = -Nhy; jh <= Nhy; jh++ ) {
int j2 = imfilter_index( j1+jh, Ny, BCy );
for ( int ih = -Nhx; ih <= Nhx; ih++ ) {
int i2 = imfilter_index( i1+ih, Nx, BCx );
bool fixed = i2 == -1 || j2 == -1 || k2 == -1;
TYPE A2 = fixed ? X : A[i2 + j2 * Nx + k2 * Nx * Ny];
tmp += H[ijkh] * A2;
ijkh++;
}
}
}
B[i1 + j1 * Nx + k1 * Nx * Ny] = tmp;
}
}
}
PROFILE_STOP( "imfilter_3D" );
}
/********************************************************
* Perform N-D filtering *
********************************************************/
template<class TYPE>
Array<TYPE> imfilter::imfilter( const Array<TYPE>& A,
const Array<TYPE>& H, const std::vector<imfilter::BC>& BC, const TYPE X )
{
IMFILTER_ASSERT( A.ndim() == H.ndim() );
IMFILTER_ASSERT( A.ndim() == BC.size() );
std::vector<size_t> Nh = H.size();
for (int d=0; d<A.ndim(); d++) {
Nh[d] = (H.size(d)-1)/2;
IMFILTER_INSIST(2*Nh[d]+1==H.size(d),"Filter must be of size 2*N+1");
}
auto B = A;
if ( A.ndim() == 1 ) {
PROFILE_START( "imfilter_1D" );
filter_direction( 1, A.size(0), 1, Nh[0], H.data(), BC[0], X, B.data() );
PROFILE_STOP( "imfilter_1D" );
} else if ( A.ndim() == 2 ) {
imfilter_2D( A.size(0), A.size(1), A.data(), Nh[0], Nh[1], H.data(), BC[0], BC[1], X, B.data() );
} else if ( A.ndim() == 3 ) {
imfilter_3D( A.size(0), A.size(1), A.size(2), A.data(),
Nh[0], Nh[1], Nh[2], H.data(), BC[0], BC[1], BC[2], X, B.data() );
} else {
IMFILTER_ERROR( "Arbitrary dimension not yet supported" );
}
return B;
}
template<class TYPE>
Array<TYPE> imfilter::imfilter( const Array<TYPE>& A, const std::vector<int>& Nh0,
std::function<TYPE(const Array<TYPE>&)> H,
const std::vector<imfilter::BC>& BC0, const TYPE X )
{
PROFILE_START( "imfilter (lambda)" );
IMFILTER_ASSERT( A.ndim() == Nh0.size() );
IMFILTER_ASSERT( A.ndim() == BC0.size() );
std::vector<size_t> Nh2( A.size() );
for (int d=0; d<A.ndim(); d++)
Nh2[d] = 2*Nh0[d]+1;
auto B = A;
Array<TYPE> data(Nh2);
IMFILTER_INSIST(A.ndim()<=3,"Not programmed for more than 3 dimensions yet");
auto N = A.size();
auto Nh = Nh0;
auto BC = BC0;
N.resize(3,1);
Nh.resize(3,0);
BC.resize(3,imfilter::BC::fixed);
for ( int k1 = 0; k1 < N[2]; k1++ ) {
for ( int j1 = 0; j1 < N[1]; j1++ ) {
for ( int i1 = 0; i1 < N[0]; i1++ ) {
for ( int kh = -Nh[2]; kh <= Nh[2]; kh++ ) {
int k2 = imfilter_index( k1+kh, N[2], BC[2] );
for ( int jh = -Nh[1]; jh <= Nh[1]; jh++ ) {
int j2 = imfilter_index( j1+jh, N[1], BC[1] );
for ( int ih = -Nh[0]; ih <= Nh[0]; ih++ ) {
int i2 = imfilter_index( i1+ih, N[0], BC[0] );
bool fixed = i2 == -1 || j2 == -1 || k2 == -1;
data(ih+Nh[0],jh+Nh[1],kh+Nh[2]) = fixed ? X : A(i2,j2,k2);
}
}
}
B(i1,j1,k1) = H( data );
}
}
}
PROFILE_STOP( "imfilter (lambda)" );
return B;
}
/********************************************************
* imfilter with separable filter functions *
********************************************************/
template<class TYPE>
Array<TYPE> imfilter::imfilter_separable( const Array<TYPE>& A,
const std::vector<Array<TYPE>>& H,
const std::vector<imfilter::BC>& boundary, const TYPE X )
{
PROFILE_START( "imfilter_separable" );
IMFILTER_ASSERT( A.ndim() == H.size() );
IMFILTER_ASSERT( A.ndim() == boundary.size() );
std::vector<size_t> Nh( H.size() );
for (int d=0; d<A.ndim(); d++) {
IMFILTER_ASSERT(H[d].ndim()==1);
Nh[d] = (H[d].length()-1)/2;
IMFILTER_INSIST(2*Nh[d]+1==H[d].length(),"Filter must be of size 2*N+1");
}
auto B = A;
for ( int d = 0; d < A.ndim(); d++ ) {
int N = A.size(d);
int Ns = 1;
int Ne = 1;
for ( int d2 = 0; d2 < d; d2++ )
Ns *= A.size(d2);
for ( int d2 = d+1; d2 < A.ndim(); d2++ )
Ne *= A.size(d2);
filter_direction( Ns, N, Ne, Nh[d], H[d].data(), boundary[d], X, B.data() );
}
PROFILE_STOP( "imfilter_separable" );
return B;
}
template<class TYPE>
Array<TYPE> imfilter::imfilter_separable( const Array<TYPE>& A, const std::vector<int>& Nh,
std::vector<std::function<TYPE(const Array<TYPE>&)>> H,
const std::vector<imfilter::BC>& boundary, const TYPE X )
{
PROFILE_START( "imfilter_separable (lambda)" );
IMFILTER_ASSERT( A.ndim() == boundary.size() );
auto B = A;
for ( int d = 0; d < A.ndim(); d++ ) {
int N = A.size(d);
int Ns = 1;
int Ne = 1;
for ( int d2 = 0; d2 < d; d2++ )
Ns *= A.size(d2);
for ( int d2 = d+1; d2 < A.ndim(); d2++ )
Ne *= A.size(d2);
filter_direction( Ns, N, Ne, Nh[d], H[d], boundary[d], X, B.data() );
}
PROFILE_STOP( "imfilter_separable (lambda)" );
return B;
}
template<class TYPE>
Array<TYPE> imfilter::imfilter_separable( const Array<TYPE>& A, const std::vector<int>& Nh,
std::vector<std::function<TYPE(int, const TYPE*)>> H,
const std::vector<imfilter::BC>& boundary, const TYPE X )
{
PROFILE_START( "imfilter_separable (function)" );
IMFILTER_ASSERT( A.ndim() == boundary.size() );
auto B = A;
for ( int d = 0; d < A.ndim(); d++ ) {
int N = A.size(d);
int Ns = 1;
int Ne = 1;
for ( int d2 = 0; d2 < d; d2++ )
Ns *= A.size(d2);
for ( int d2 = d+1; d2 < A.ndim(); d2++ )
Ne *= A.size(d2);
filter_direction( Ns, N, Ne, Nh[d], H[d], boundary[d], X, B.data() );
}
PROFILE_STOP( "imfilter_separable (function)" );
return B;
}

4
common/pmmc.h
View File

@@ -6,9 +6,9 @@
#include <iostream>
#include <fstream>
#include <math.h>
#include "Array.h"
#include "common/Array.h"
#include "PointList.h"
#include "Utilities.h"
#include "common/Utilities.h"
using namespace std;

6
tests/BlobAnalysis.cpp
View File

@@ -436,7 +436,7 @@ int main(int argc, char **argv)
*/
FILE *PHASE = fopen("Phase.dat","wb");
fwrite(Phase.get(),8,Nx*Ny*Nz,PHASE);
fwrite(Phase.data(),8,Nx*Ny*Nz,PHASE);
fclose(PHASE);
// Compute the porosity
@@ -998,12 +998,12 @@ int main(int argc, char **argv)
FILE *BLOBS;
BLOBS = fopen("Blobs.dat","wb");
fwrite(LocalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
fwrite(LocalBlobID.data(),4,Nx*Ny*Nz,BLOBS);
fclose(BLOBS);
FILE *DISTANCE;
DISTANCE = fopen("SignDist.dat","wb");
fwrite(SignDist.get(),8,Nx*Ny*Nz,DISTANCE);
fwrite(SignDist.data(),8,Nx*Ny*Nz,DISTANCE);
fclose(DISTANCE);
}

6
tests/BlobIdentify.cpp
View File

@@ -141,8 +141,8 @@ void readRankData( int proc, int nx, int ny, int nz, DoubleArray& Phase, DoubleA
char file1[40], file2[40];
sprintf(file1,"SignDist.%05d",proc);
sprintf(file2,"Phase.%05d",proc);
ReadBinaryFile(file1, Phase.get(), nx*ny*nz);
ReadBinaryFile(file2, SignDist.get(), nx*ny*nz);
ReadBinaryFile(file1, Phase.data(), nx*ny*nz);
ReadBinaryFile(file2, SignDist.data(), nx*ny*nz);
}
@@ -316,7 +316,7 @@ int main(int argc, char **argv)
FILE *BLOBS = fopen("Blobs.dat","wb");
fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
fwrite(GlobalBlobID.data(),4,Nx*Ny*Nz,BLOBS);
fclose(BLOBS);
MPI_Finalize();
return 0;

8
tests/BlobIdentifyParallel.cpp
View File

@@ -39,8 +39,8 @@ void readRankData( int proc, int nx, int ny, int nz, DoubleArray& Phase, DoubleA
char file1[40], file2[40];
sprintf(file1,"SignDist.%05d",proc);
sprintf(file2,"Phase.%05d",proc);
ReadBinaryFile(file1, Phase.get(), nx*ny*nz);
ReadBinaryFile(file2, SignDist.get(), nx*ny*nz);
ReadBinaryFile(file1, Phase.data(), nx*ny*nz);
ReadBinaryFile(file2, SignDist.data(), nx*ny*nz);
}
@@ -117,13 +117,13 @@ int main(int argc, char **argv)
char LocalRankFilename[100];
sprintf(LocalRankFilename,"BlobLabel.%05i",rank);
FILE *BLOBLOCAL = fopen(LocalRankFilename,"wb");
fwrite(GlobalBlobID.get(),4,GlobalBlobID.length(),BLOBLOCAL);
fwrite(GlobalBlobID.data(),4,GlobalBlobID.length(),BLOBLOCAL);
fclose(BLOBLOCAL);
printf("Wrote BlobLabel.%05i \n",rank);
/*FILE *BLOBS = fopen("Blobs.dat","wb");
fwrite(GlobalBlobID.get(),4,Nx*Ny*Nz,BLOBS);
fwrite(GlobalBlobID.data(),4,Nx*Ny*Nz,BLOBS);
fclose(BLOBS);*/
#ifdef PROFILE
PROFILE_STOP("main");

4
tests/ColorToBinary.cpp
View File

@@ -266,11 +266,11 @@ int main(int argc, char **argv)
fclose(OUTFILE);
OUTFILE = fopen("Phase.dat","wb");
fwrite(Phase.get(),8,Nx*Ny*Nz,OUTFILE);
fwrite(Phase.data(),8,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
/* OUTFILE = fopen("SignDist.dat","wb");
fwrite(SignDist.get(),8,Nx*Ny*Nz,OUTFILE);
fwrite(SignDist.data(),8,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
*/

6
tests/TestBlobAnalyze.cpp
View File

@@ -20,8 +20,8 @@ void readRankData( int proc, int nx, int ny, int nz, DoubleArray& Phase, DoubleA
char file1[40], file2[40];
sprintf(file1,"SignDist.%05d",proc);
sprintf(file2,"Phase.%05d",proc);
ReadBinaryFile(file1, Phase.get(), nx*ny*nz);
ReadBinaryFile(file2, SignDist.get(), nx*ny*nz);
ReadBinaryFile(file1, Phase.data(), nx*ny*nz);
ReadBinaryFile(file2, SignDist.data(), nx*ny*nz);
}
inline void WriteBlobs(TwoPhase Averages){
printf("Writing the blob list \n");
@@ -262,7 +262,7 @@ int main(int argc, char **argv)
//...........................................................................
MPI_Barrier(comm);
//.......................................................................
SignedDistance(Averages.Phase.get(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
SignedDistance(Averages.Phase.data(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
Dm.iproc,Dm.jproc,Dm.kproc,Dm.nprocx,Dm.nprocy,Dm.nprocz);
//.......................................................................
// Assign the phase ID field based on the signed distance

28
tests/TestBubble.cpp
View File

@@ -1198,8 +1198,8 @@ int main(int argc, char **argv)
//...........................................................................
InitD3Q19(ID, f_even, f_odd, Nx, Ny, Nz);
//......................................................................
// InitDenColorDistance(ID, Copy, Phi, SDs.get(), das, dbs, beta, xIntPos, Nx, Ny, Nz, S);
InitDenColorDistance(ID, Den, Phi, SDs.get(), das, dbs, beta, xIntPos, Nx, Ny, Nz);
// InitDenColorDistance(ID, Copy, Phi, SDs.data(), das, dbs, beta, xIntPos, Nx, Ny, Nz, S);
InitDenColorDistance(ID, Den, Phi, SDs.data(), das, dbs, beta, xIntPos, Nx, Ny, Nz);
InitD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
//......................................................................
@@ -1212,7 +1212,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
WriteLocalSolidID(LocalRankFilename, id, N);
sprintf(LocalRankFilename,"%s%s","SDs.",LocalRankString);
WriteLocalSolidDistance(LocalRankFilename, SDs.get(), N);
WriteLocalSolidDistance(LocalRankFilename, SDs.data(), N);
//.......................................................................
if (Restart == true){
if (rank==0) printf("Reading restart file! \n");
@@ -1345,7 +1345,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Phase_tplus.get(),Phi,N*sizeof(double));
CopyToHost(Phase_tplus.data(),Phi,N*sizeof(double));
//...........................................................................
//...........................................................................
// Copy the data for for the analysis timestep
@@ -1354,8 +1354,8 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Phase.get(),Phi,N*sizeof(double));
CopyToHost(Press.get(),Pressure,N*sizeof(double));
CopyToHost(Phase.data(),Phi,N*sizeof(double));
CopyToHost(Press.data(),Pressure,N*sizeof(double));
MPI_Barrier(comm);
//...........................................................................
@@ -1774,14 +1774,14 @@ int main(int argc, char **argv)
// Copy the phase indicator field for the later timestep
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Phase_tminus.get(),Phi,N*sizeof(double));
CopyToHost(Phase_tplus.get(),Phi,N*sizeof(double));
CopyToHost(Phase.get(),Phi,N*sizeof(double));
CopyToHost(Press.get(),Pressure,N*sizeof(double));
CopyToHost(Phase_tminus.data(),Phi,N*sizeof(double));
CopyToHost(Phase_tplus.data(),Phi,N*sizeof(double));
CopyToHost(Phase.data(),Phi,N*sizeof(double));
CopyToHost(Press.data(),Pressure,N*sizeof(double));
double temp=0.5/beta;
for (n=0; n<N; n++){
double value = Phase.get()[n];
double value = Phase.data()[n];
SDn(n) = temp*log((1.0+value)/(1.0-value));
}
@@ -2218,12 +2218,12 @@ int main(int argc, char **argv)
//************************************************************************/
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
// printf("Local File Name = %s \n",LocalRankFilename);
// CopyToHost(Phase.get(),Phi,N*sizeof(double));
// CopyToHost(Phase.data(),Phi,N*sizeof(double));
FILE *PHASE;
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Press.get(),8,N,PHASE);
// fwrite(MeanCurvature.get(),8,N,PHASE);
fwrite(Press.data(),8,N,PHASE);
// fwrite(MeanCurvature.data(),8,N,PHASE);
fclose(PHASE);
/* double *DensityValues;

4
tests/TestSphereCurvature.cpp
View File

@@ -79,7 +79,7 @@ int main(int argc, char **argv)
DoubleArray CubeValues(2,2,2);
// Compute the signed distance function
SignedDistance(Phase.get(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
SignedDistance(Phase.data(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
@@ -89,7 +89,7 @@ int main(int argc, char **argv)
}
}
}
SignedDistance(SignDist.get(),0,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
SignedDistance(SignDist.data(),0,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
pmmc_MeshCurvature(Phase, MeanCurvature, GaussCurvature, Nx, Ny, Nz);

2
tests/TestTwoPhase.cpp
View File

@@ -106,7 +106,7 @@ int main(int argc, char **argv)
if (rank==0){
FILE *PHASE = fopen("Phase.00000","wb");
fwrite(Averages.MeanCurvature.get(),8,Nx*Ny*Nz,PHASE);
fwrite(Averages.MeanCurvature.data(),8,Nx*Ny*Nz,PHASE);
fclose(PHASE);
}
// ****************************************************

8
tests/lbpm_BlobAnalysis.cpp
View File

@@ -121,16 +121,16 @@ int main(int argc, char **argv)
if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
//.......................................................................
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Averages.SDs.get(), N);
ReadBinaryFile(LocalRankFilename, Averages.SDs.data(), N);
MPI_Barrier(comm);
// sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
//ReadBinaryFile(LocalRankFilename, Averages.Press.get(), N);
//ReadBinaryFile(LocalRankFilename, Averages.Press.data(), N);
//MPI_Barrier(comm);
if (rank == 0) cout << "Domain set." << endl;
//.......................................................................
sprintf(LocalRankFilename,"%s%s","Label_NWP.",LocalRankString);
ReadBlobFile(LocalRankFilename, Averages.Label_NWP.get(), N);
ReadBlobFile(LocalRankFilename, Averages.Label_NWP.data(), N);
MPI_Barrier(comm);
if (rank == 0) cout << "Label_NWP set." << endl;
@@ -227,7 +227,7 @@ int main(int argc, char **argv)
}
// BlobContainer Blobs;
DoubleArray RecvBuffer(dimx);
// MPI_Allreduce(&Averages.ComponentAverages_NWP.get(),&Blobs.get(),1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
// MPI_Allreduce(&Averages.ComponentAverages_NWP.data(),&Blobs.data(),1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Barrier(comm);
if (rank==0) printf("All ranks passed gate \n");

2
tests/lbpm_captube_pp.cpp
View File

@@ -198,7 +198,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.get(),8,Averages.SDs.length(),DIST);
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fclose(DIST);
sprintf(LocalRankFilename,"ID.%05i",rank);

24
tests/lbpm_color_simulator.cpp
View File

@@ -372,7 +372,7 @@ int main(int argc, char **argv)
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// WriteLocalSolidID(LocalRankFilename, id, N);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Averages->SDs.get(), N);
ReadBinaryFile(LocalRankFilename, Averages->SDs.data(), N);
MPI_Barrier(comm);
if (rank == 0) cout << "Domain set." << endl;
@@ -553,7 +553,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy signed distance for device initialization
CopyToDevice(dvcSignDist, Averages->SDs.get(), dist_mem_size);
CopyToDevice(dvcSignDist, Averages->SDs.data(), dist_mem_size);
//...........................................................................
int logcount = 0; // number of surface write-outs
@@ -686,7 +686,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy the phase indicator field for the earlier timestep
DeviceBarrier();
CopyToHost(Averages->Phase_tplus.get(),Phi,N*sizeof(double));
CopyToHost(Averages->Phase_tplus.data(),Phi,N*sizeof(double));
//...........................................................................
//...........................................................................
// Copy the data for for the analysis timestep
@@ -695,11 +695,11 @@ int main(int argc, char **argv)
//...........................................................................
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
CopyToHost(Averages->Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages->Press.get(),Pressure,N*sizeof(double));
CopyToHost(Averages->Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages->Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages->Vel_z.get(),&Velocity[2*N],N*sizeof(double));
CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
CopyToHost(Averages->Press.data(),Pressure,N*sizeof(double));
CopyToHost(Averages->Vel_x.data(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages->Vel_y.data(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages->Vel_z.data(),&Velocity[2*N],N*sizeof(double));
//...........................................................................
if (rank==0) printf("********************************************************\n");
@@ -934,7 +934,7 @@ int main(int argc, char **argv)
int NumberComponents_NWP = ComputeGlobalPhaseComponent(Mask.Nx-2,Mask.Ny-2,Mask.Nz-2,Mask.rank_info,Averages->PhaseID,1,Averages->Label_NWP);
printf("Number of non-wetting phase components: %i \n ",NumberComponents_NWP);
DeviceBarrier();
CopyToHost(Averages->Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
*/
/* Averages->WriteSurfaces(0);
@@ -942,17 +942,17 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
FILE *PHASE;
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Averages->SDn.get(),8,N,PHASE);
fwrite(Averages->SDn.data(),8,N,PHASE);
fclose(PHASE);
*/
/* sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
FILE *PRESS;
PRESS = fopen(LocalRankFilename,"wb");
fwrite(Averages->Press.get(),8,N,PRESS);
fwrite(Averages->Press.data(),8,N,PRESS);
fclose(PRESS);
CopyToHost(Averages->Phase.get(),Phi,N*sizeof(double));
CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
double * Grad;
Grad = new double [3*N];
CopyToHost(Grad,ColorGrad,3*N*sizeof(double));

23
tests/lbpm_color_simulator.h
View File

@@ -12,6 +12,13 @@ enum AnalysisType{ AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02,
CopySimState=0x04, ComputeAverages=0x08, CreateRestart=0x10, WriteVis=0x20 };
template<class TYPE>
void DeleteArray( const TYPE *p )
{
delete [] p;
}
// Structure used to store ids
struct AnalysisWaitIdStruct {
ThreadPool::thread_id_t blobID;
@@ -279,14 +286,14 @@ void run_analysis( int timestep, int restart_interval,
(type&CopySimState)!=0 || (type&IdentifyBlobs)!=0 )
{
phase = std::shared_ptr<DoubleArray>(new DoubleArray(Nx,Ny,Nz));
CopyToHost(phase->get(),Phi,N*sizeof(double));
CopyToHost(phase->data(),Phi,N*sizeof(double));
}
if ( (type&CopyPhaseIndicator)!=0 ) {
memcpy(Averages.Phase_tplus.get(),phase->get(),N*sizeof(double));
memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double));
//Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
}
if ( (type&ComputeAverages)!=0 ) {
memcpy(Averages.Phase_tminus.get(),phase->get(),N*sizeof(double));
memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double));
//Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tminus);
}
if ( (type&CopySimState) != 0 ) {
@@ -301,11 +308,11 @@ void run_analysis( int timestep, int restart_interval,
tpool.wait(wait.vis); // Make sure we are done using analysis before modifying
PROFILE_STOP("Copy-Wait",1);
PROFILE_START("Copy-State",1);
memcpy(Averages.Phase.get(),phase->get(),N*sizeof(double));
CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
CopyToHost(Averages.Press.data(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.data(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.data(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.data(),&Velocity[2*N],N*sizeof(double));
PROFILE_STOP("Copy-State",1);
}
std::shared_ptr<double> cDen, cDistEven, cDistOdd;

4
tests/lbpm_disc_pp.cpp
View File

@@ -294,7 +294,7 @@ int main(int argc, char **argv)
if (nprocx > 1 && rank==0) printf("Disc packs are 2D -- are you sure you want nprocx > 1? \n");
//.......................................................................
SignedDistanceDiscPack(SignDist.get(),ndiscs,cx,cy,rad,Lx,Ly,Lz,Nx,Ny,Nz,
SignedDistanceDiscPack(SignDist.data(),ndiscs,cx,cy,rad,Lx,Ly,Lz,Nx,Ny,Nz,
iproc,jproc,kproc,nprocx,nprocy,nprocz);
//.......................................................................
// Assign walls in the signed distance functions (x,y boundaries)
@@ -375,7 +375,7 @@ int main(int argc, char **argv)
//.......................................................................
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
WriteLocalSolidDistance(LocalRankFilename, SignDist.get(), N);
WriteLocalSolidDistance(LocalRankFilename, SignDist.data(), N);
//......................................................................
// ****************************************************

2
tests/lbpm_morphdrain_pp.cpp
View File

@@ -133,7 +133,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"rb");
size_t ReadSignDist;
ReadSignDist=fread(SignDist.get(),8,N,DIST);
ReadSignDist=fread(SignDist.data(),8,N,DIST);
if (ReadSignDist != size_t(N)) printf("lbpm_morphdrain_pp: Error reading signed distance function (rank=%i)\n",rank);
fclose(DIST);

12
tests/lbpm_permeability_simulator.cpp
View File

@@ -282,7 +282,7 @@ int main(int argc, char **argv)
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// WriteLocalSolidID(LocalRankFilename, id, N);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Averages.SDs.get(), N);
ReadBinaryFile(LocalRankFilename, Averages.SDs.data(), N);
MPI_Barrier(comm);
if (rank == 0) cout << "Domain set." << endl;
@@ -402,7 +402,7 @@ int main(int argc, char **argv)
//...........................................................................
// Copy signed distance for device initialization
CopyToDevice(dvcSignDist, Averages.SDs.get(), dist_mem_size);
CopyToDevice(dvcSignDist, Averages.SDs.data(), dist_mem_size);
//...........................................................................
int logcount = 0; // number of surface write-outs
@@ -509,10 +509,10 @@ int main(int argc, char **argv)
DeviceBarrier();
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
ComputeVelocityD3Q19(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
CopyToHost(Averages.Press.data(),Pressure,N*sizeof(double));
CopyToHost(Averages.Vel_x.data(),&Velocity[0],N*sizeof(double));
CopyToHost(Averages.Vel_y.data(),&Velocity[N],N*sizeof(double));
CopyToHost(Averages.Vel_z.data(),&Velocity[2*N],N*sizeof(double));
// Way more work than necessary -- this is just to get the solid interfacial area!!
Averages.Initialize();

2
tests/lbpm_random_pp.cpp
View File

@@ -148,7 +148,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"rb");
size_t ReadSignDist;
ReadSignDist=fread(SignDist.get(),8,N,DIST);
ReadSignDist=fread(SignDist.data(),8,N,DIST);
if (ReadSignDist != size_t(N)) printf("lbpm_random_pp: Error reading signed distance function (rank=%i)\n",rank);
fclose(DIST);

2
tests/lbpm_segmented_pp.cpp
View File

@@ -332,7 +332,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.get(),8,Averages.SDs.length(),DIST);
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fclose(DIST);
/* // Solve for the position of the non-wetting phase

4
tests/lbpm_sphere_pp.cpp
View File

@@ -185,7 +185,7 @@ int main(int argc, char **argv)
MPI_Bcast(&D,1,MPI_DOUBLE,0,comm);
//.......................................................................
SignedDistance(SignDist.get(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
SignedDistance(SignDist.data(),nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,
iproc,jproc,kproc,nprocx,nprocy,nprocz);
//.......................................................................
// Assign the phase ID field based on the signed distance
@@ -242,7 +242,7 @@ int main(int argc, char **argv)
//.......................................................................
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
WriteLocalSolidDistance(LocalRankFilename, SignDist.get(), N);
WriteLocalSolidDistance(LocalRankFilename, SignDist.data(), N);
//......................................................................
// ****************************************************

2
tests/lbpm_squaretube_pp.cpp
View File

@@ -198,7 +198,7 @@ int main(int argc, char **argv)
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"wb");
fwrite(Averages.SDs.get(),8,Averages.SDs.length(),DIST);
fwrite(Averages.SDs.data(),8,Averages.SDs.length(),DIST);
fclose(DIST);
sprintf(LocalRankFilename,"ID.%05i",rank);

1046
tests/lbpm_uCT_pp.cpp
View File

File diff suppressed because it is too large Load Diff