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

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

Browse Source
This commit is contained in:
James McClure
2019-12-03 13:58:09 -05:00
parent cf5a284f6d 442d42bb91
commit 9ea68b848c
14 changed files with 458 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files

52
analysis/SubPhase.cpp
View File

@@ -109,35 +109,35 @@ SubPhase::~SubPhase()
void SubPhase::Write(int timestep)
{
if (Dm->rank()==0){
fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.p, gwd.p, gnc.p, gnd.p);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.V, gwc.A, gwc.H, gwc.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gnc.V, gnc.A, gnc.H, gnc.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gnd.V, gnd.A, gnd.H, gnd.X, gnd.Nc);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",giwn.V, giwn.A, giwn.H, giwn.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i\n",giwnc.V, giwnc.A, giwnc.H, giwnc.X, giwnc.Nc);
fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.p, gwd.p, gnc.p, gnd.p);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.V, gwc.A, gwc.H, gwc.X);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gnc.V, gnc.A, gnc.H, gnc.X);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",gnd.V, gnd.A, gnd.H, gnd.X, gnd.Nc);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",giwn.V, giwn.A, giwn.H, giwn.X);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i\n",giwnc.V, giwnc.A, giwnc.H, giwnc.X, giwnc.Nc);
fflush(SUBPHASE);
}
else{
fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",wc.p, wd.p, nc.p, nd.p);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",wc.V, wc.A, wc.H, wc.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",nc.V, nc.A, nc.H, nc.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",iwn.V, iwn.A, iwn.H, iwn.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g\n",iwnc.V, iwnc.A, iwnc.H, iwnc.X);
fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.p, wd.p, nc.p, nd.p);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.V, wc.A, wc.H, wc.X);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",nc.V, nc.A, nc.H, nc.X);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",iwn.V, iwn.A, iwn.H, iwn.X);
fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g\n",iwnc.V, iwnc.A, iwnc.H, iwnc.X);
}
}

8
common/Array.hpp
View File

@@ -694,15 +694,15 @@ Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::repmat(
for ( size_t i4 = 0, index = 0; i4 < N1[4]; i4++ ) {
for ( size_t j4 = 0; j4 < Nr[4]; j4++ ) {
for ( size_t i3 = 0; i3 < N1[3]; i3++ ) {
for ( size_t j4 = 0; j4 < Nr[3]; j4++ ) {
for ( size_t j3 = 0; j3 < Nr[3]; j3++ ) {
for ( size_t i2 = 0; i2 < N1[2]; i2++ ) {
for ( size_t j4 = 0; j4 < Nr[2]; j4++ ) {
for ( size_t j2 = 0; j2 < Nr[2]; j2++ ) {
for ( size_t i1 = 0; i1 < N1[1]; i1++ ) {
for ( size_t j4 = 0; j4 < Nr[1]; j4++ ) {
for ( size_t j1 = 0; j1 < Nr[1]; j1++ ) {
for ( size_t i0 = 0; i0 < N1[0]; i0++ ) {
size_t k = d_size.index( i0, i1, i2, i3, i4 );
TYPE x = d_data[k];
for ( size_t j4 = 0; j4 < Nr[0]; j4++, index++ )
for ( size_t j0 = 0; j0 < Nr[0]; j0++, index++ )
y2[index] = x;
}
}

53
common/Communication.cpp
View File

@@ -4,16 +4,28 @@
/********************************************************
* Structure to store the rank info *
********************************************************/
inline int getRankForBlock( int nprocx, int nprocy, int nprocz, int i, int j, int k )
int RankInfoStruct::getRankForBlock( int i, int j, int k ) const
{
int i2 = (i+nprocx)%nprocx;
int j2 = (j+nprocy)%nprocy;
int k2 = (k+nprocz)%nprocz;
return i2 + j2*nprocx + k2*nprocx*nprocy;
int i2 = (i+nx)%nx;
int j2 = (j+ny)%ny;
int k2 = (k+nz)%nz;
return i2 + j2*nx + k2*nx*ny;
}
RankInfoStruct::RankInfoStruct()
{
memset(this,0,sizeof(RankInfoStruct));
nx = 0;
ny = 0;
nz = 0;
ix = -1;
jy = -1;
kz = -1;
for (int i=-1; i<=1; i++) {
for (int j=-1; j<=1; j++) {
for (int k=-1; k<=1; k++) {
rank[i+1][j+1][k+1] = -1;
}
}
}
}
RankInfoStruct::RankInfoStruct( int rank0, int nprocx, int nprocy, int nprocz )
{
@@ -21,17 +33,30 @@ RankInfoStruct::RankInfoStruct( int rank0, int nprocx, int nprocy, int nprocz )
nx = nprocx;
ny = nprocy;
nz = nprocz;
ix = rank0%nprocx;
jy = (rank0/nprocx)%nprocy;
kz = rank0/(nprocx*nprocy);
for (int i=-1; i<=1; i++) {
for (int j=-1; j<=1; j++) {
for (int k=-1; k<=1; k++) {
rank[i+1][j+1][k+1] = getRankForBlock(nprocx,nprocy,nprocz,ix+i,jy+j,kz+k);
if ( rank0 >= nprocx * nprocy * nprocz ) {
ix = -1;
jy = -1;
kz = -1;
for (int i=-1; i<=1; i++) {
for (int j=-1; j<=1; j++) {
for (int k=-1; k<=1; k++) {
rank[i+1][j+1][k+1] = -1;
}
}
}
} else {
ix = rank0%nprocx;
jy = (rank0/nprocx)%nprocy;
kz = rank0/(nprocx*nprocy);
for (int i=-1; i<=1; i++) {
for (int j=-1; j<=1; j++) {
for (int k=-1; k<=1; k++) {
rank[i+1][j+1][k+1] = getRankForBlock(ix+i,jy+j,kz+k);
}
}
}
ASSERT(rank[1][1][1]==rank0);
}
ASSERT(rank[1][1][1]==rank0);
}

7
common/Communication.h
View File

@@ -31,9 +31,16 @@ struct RankInfoStruct {
int rank[3][3][3]; //!< The rank for the neighbor [i][j][k]
RankInfoStruct();
RankInfoStruct( int rank, int nprocx, int nprocy, int nprocz );
int getRankForBlock( int i, int j, int k ) const;
};
//! Redistribute domain data (dst may be smaller than the src)
template<class TYPE>
Array<TYPE> redistribute( const RankInfoStruct& src_rank, const Array<TYPE>& src_data,
const RankInfoStruct& dst_rank, std::array<int,3> dst_size, MPI_Comm comm );
/*!
* @brief Communicate halo
* @details Fill the halo cells in an array from the neighboring processes

85
common/Communication.hpp
View File

@@ -8,7 +8,90 @@
/********************************************************
* Structure to store the rank info *
* Redistribute data between two grids *
********************************************************/
template<class TYPE>
Array<TYPE> redistribute( const RankInfoStruct& src_rank, const Array<TYPE>& src_data,
const RankInfoStruct& dst_rank, std::array<int,3> dst_size, MPI_Comm comm )
{
#ifdef USE_MPI
// Get the src size
std::array<int,3> src_size;
int size0[3] = { (int) src_data.size(0), (int) src_data.size(1), (int) src_data.size(2) };
MPI_Allreduce( size0, src_size.data(), 3, MPI_INT, MPI_MAX, comm );
if ( !src_data.empty() )
ASSERT( src_size[0] == size0[0] && src_size[1] == size0[1] && src_size[2] == size0[2] );
// Check that dst_size matches on all ranks
MPI_Allreduce( dst_size.data(), size0, 3, MPI_INT, MPI_MAX, comm );
ASSERT( dst_size[0] == size0[0] && dst_size[1] == size0[1] && dst_size[2] == size0[2] );
// Function to get overlap range
auto calcOverlap = []( int i1[3], int i2[3], int j1[3], int j2[3] ) {
std::vector<size_t> index;
if ( i1[0] > j2[0] || i2[0] < j1[0] || i1[1] > j2[1] || i2[1] < j1[1] || i1[2] > j2[2] || i2[2] < j1[2] )
return index;
index.resize( 6 );
index[0] = std::max( j1[0] - i1[0], 0 );
index[1] = std::min( j2[0] - i1[0], i2[0] - i1[0] );
index[2] = std::max( j1[1] - i1[1], 0 );
index[3] = std::min( j2[1] - i1[1], i2[1] - i1[1] );
index[4] = std::max( j1[2] - i1[2], 0 );
index[5] = std::min( j2[2] - i1[2], i2[2] - i1[2] );
return index;
};
// Pack and send my data to the appropriate ranks (including myself)
std::vector<int> send_rank;
std::vector<Array<TYPE>> send_data;
if ( !src_data.empty() ) {
int i1[3] = { src_size[0] * src_rank.ix, src_size[1] * src_rank.jy, src_size[2] * src_rank.kz };
int i2[3] = { i1[0] + src_size[0] - 1, i1[1] + src_size[1] - 1, i1[2] + src_size[2] - 1 };
for ( size_t i=0; i<dst_rank.nx; i++ ) {
for ( size_t j=0; j<dst_rank.ny; j++ ) {
for ( size_t k=0; k<dst_rank.nz; k++ ) {
int j1[3] = { i * dst_size[0], j * dst_size[1], k * dst_size[2] };
int j2[3] = { j1[0] + dst_size[0] - 1, j1[1] + dst_size[1] - 1, j1[2] + dst_size[2] - 1 };
auto index = calcOverlap( i1, i2, j1, j2 );
if ( index.empty() )
continue;
send_rank.push_back( dst_rank.getRankForBlock(i,j,k) );
send_data.push_back( src_data.subset( index ) );
}
}
}
}
std::vector<MPI_Request> send_request( send_rank.size() );
for (size_t i=0; i<send_rank.size(); i++)
MPI_Isend( send_data[i].data(), sizeof(TYPE)*send_data[i].length(), MPI_BYTE, send_rank[i], 5462, comm, &send_request[i]);
// Unpack data from the appropriate ranks (including myself)
Array<TYPE> dst_data( dst_size[0], dst_size[1], dst_size[2] );
int i1[3] = { dst_size[0] * dst_rank.ix, dst_size[1] * dst_rank.jy, dst_size[2] * dst_rank.kz };
int i2[3] = { i1[0] + dst_size[0] - 1, i1[1] + dst_size[1] - 1, i1[2] + dst_size[2] - 1 };
for ( size_t i=0; i<src_rank.nx; i++ ) {
for ( size_t j=0; j<src_rank.ny; j++ ) {
for ( size_t k=0; k<src_rank.nz; k++ ) {
int j1[3] = { i * src_size[0], j * src_size[1], k * src_size[2] };
int j2[3] = { j1[0] + src_size[0] - 1, j1[1] + src_size[1] - 1, j1[2] + src_size[2] - 1 };
auto index = calcOverlap( i1, i2, j1, j2 );
if ( index.empty() )
continue;
int rank = src_rank.getRankForBlock(i,j,k);
Array<TYPE> data( index[1] - index[0] + 1, index[3] - index[2] + 1, index[5] - index[4] + 1 );
MPI_Recv( data.data(), sizeof(TYPE)*data.length(), MPI_BYTE, rank, 5462, comm, MPI_STATUS_IGNORE );
dst_data.copySubset( index, data );
}
}
}
// Free data
MPI_Waitall( send_request.size(), send_request.data(), MPI_STATUSES_IGNORE );
return dst_data;
#else
return src_data.subset( { 0, dst_size[0]-1, 0, dst_size[1]-1, 0, dst_size[2]-1 );
#endif
}
/********************************************************
* Structure to fill halo cells *
********************************************************/
template<class TYPE>
fillHalo<TYPE>::fillHalo( MPI_Comm comm_, const RankInfoStruct& info_,

102
common/ReadMicroCT.cpp Normal file
View File

@@ -0,0 +1,102 @@
#include "common/ReadMicroCT.h"
#include "common/Utilities.h"
#include "zlib.h"
#include <cstring>
// Read a file into memory
std::vector<char> readFile( const std::string& filename )
{
auto fid = fopen( filename.c_str(), "rb" );
INSIST( fid, "File does not exist: " + filename );
fseek( fid, 0, SEEK_END );
size_t bytes = ftell(fid);
fseek( fid, 0, SEEK_SET );
std::vector<char> data( bytes );
size_t bytes2 = fread( data.data(), 1, bytes, fid );
ASSERT( bytes == bytes2 );
fclose( fid );
return data;
}
// Decompress a gzip buffer
std::vector<char> gunzip( const std::vector<char>& in )
{
z_stream stream;
std::vector<char> out( 1000000 );
stream.next_in = (Bytef*) in.data();
stream.avail_in = in.size();
stream.total_in = 0;
stream.zalloc = Z_NULL;
stream.zfree = Z_NULL;
stream.opaque = Z_NULL;
stream.next_out = (Bytef*) out.data();
stream.avail_out = out.size();
stream.total_out = 0;
ASSERT( inflateInit2(&stream, 16+MAX_WBITS) == Z_OK );
bool finished = inflate(&stream, Z_SYNC_FLUSH) == Z_STREAM_END;
while ( !finished && stream.msg == Z_NULL ) {
out.resize( 2 * out.size() );
stream.next_out = (Bytef*) &out[stream.total_out];
stream.avail_out = out.size() - stream.total_out;
finished = inflate(&stream, Z_SYNC_FLUSH) == Z_STREAM_END;
}
ASSERT( stream.msg == Z_NULL );
out.resize( stream.total_out );
inflateEnd(&stream);
return out;
}
// Read the compressed micro CT data
Array<uint8_t> readMicroCT( const std::string& filename )
{
auto in = readFile( filename );
auto out = gunzip( in );
ASSERT( out.size() == 1024*1024*1024 );
Array<uint8_t> data( 1024, 1024, 1024 );
memcpy( data.data(), out.data(), data.length() );
return data;
}
// Read the compressed micro CT data and distribute
Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm )
{
// Get the local problem info
auto n = domain.getVector<int>( "n" );
int rank = comm_rank(MPI_COMM_WORLD);
auto nproc = domain.getVector<int>( "nproc" );
RankInfoStruct rankInfo( rank, nproc[0], nproc[1], nproc[2] );
// Determine the largest file number to get
int Nfx = ( n[0] * rankInfo.nx + 1023 ) / 1024;
int Nfy = ( n[1] * rankInfo.ny + 1023 ) / 1024;
int Nfz = ( n[2] * rankInfo.nz + 1023 ) / 1024;
// Load one of the files if rank < largest file
Array<uint8_t> data;
RankInfoStruct srcRankInfo( rank, Nfx, Nfy, Nfz );
if ( srcRankInfo.ix >= 0 ) {
auto filename = domain.getScalar<std::string>( "Filename" );
char tmp[100];
if ( filename.find( "0x_0y_0z.gbd.gz" ) != std::string::npos ) {
sprintf( tmp, "%ix_%iy_%iz.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
filename = filename.replace( filename.find( "0x_0y_0z.gbd.gz" ), 15, std::string( tmp ) );
} else if ( filename.find( "x0_y0_z0.gbd.gz" ) != std::string::npos ) {
sprintf( tmp, "x%i_y%i_z%i.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
filename = filename.replace( filename.find( "x0_y0_z0.gbd.gz" ), 15, std::string( tmp ) );
} else {
ERROR( "Invalid name for first file" );
}
data = readMicroCT( filename );
}
// Redistribute the data
data = redistribute( srcRankInfo, data, rankInfo, { n[0], n[1], n[2] }, comm );
return data;
}

15
common/ReadMicroCT.h Normal file
View File

@@ -0,0 +1,15 @@
#ifndef READMICROCT_H
#define READMICROCT_H
#include "common/Array.h"
#include "common/Communication.h"
#include "common/Database.h"
Array<uint8_t> readMicroCT( const std::string& filename );
Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm );
#endif

20
cpu/Color.cpp
View File

@@ -1364,11 +1364,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, do
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
if (C==0.0) C=1.0;
nx = nx/C;
ny = ny/C;
nz = nz/C;
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=0
fq = dist[n];
rho = fq;
@@ -1949,10 +1950,11 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *di
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
if (C==0.0) C=1.0;
nx = nx/C;
ny = ny/C;
nz = nz/C;
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=0
fq = dist[n];

18
gpu/Color.cu
View File

@@ -1356,10 +1356,11 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
if (C==0.0) C=1.0;
nx = nx/C;
ny = ny/C;
nz = nz/C;
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=0
fq = dist[n];
@@ -1942,10 +1943,11 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
if (C==0.0) C=1.0;
nx = nx/C;
ny = ny/C;
nz = nz/C;
double ColorMag = C;
if (C==0.0) ColorMag=1.0;
nx = nx/ColorMag;
ny = ny/ColorMag;
nz = nz/ColorMag;
// q=0
fq = dist[n];

90
models/ColorModel.cpp
View File

@@ -508,6 +508,8 @@ void ScaLBL_ColorModel::Run(){
double NOISE_THRESHOLD = 0.0;
double BUMP_RATE = 2.0;
bool USE_BUMP_RATE = false;
int RESCALE_FORCE_COUNT = 0;
int RESCALE_FORCE_MAX = 0;
auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
if (protocol == "image sequence"){
@@ -549,6 +551,10 @@ void ScaLBL_ColorModel::Run(){
if (color_db->keyExists( "capillary_number" )){
capillary_number = color_db->getScalar<double>( "capillary_number" );
SET_CAPILLARY_NUMBER=true;
//RESCALE_FORCE_MAX = 1;
}
if (analysis_db->keyExists( "rescale_force_count" )){
RESCALE_FORCE_MAX = analysis_db->getScalar<int>( "rescale_force_count" );
}
if (color_db->keyExists( "timestep" )){
timestep = color_db->getScalar<int>( "timestep" );
@@ -591,6 +597,7 @@ void ScaLBL_ColorModel::Run(){
MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
}
if (rank==0){
printf("********************************************************\n");
if (protocol == "image sequence"){
@@ -768,6 +775,23 @@ void ScaLBL_ColorModel::Run(){
isSteady = true;
if (CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS)
isSteady = true;
if (SET_CAPILLARY_NUMBER && RESCALE_FORCE_COUNT < RESCALE_FORCE_MAX){
RESCALE_FORCE_COUNT++;
Fx *= capillary_number / Ca;
Fy *= capillary_number / Ca;
Fz *= capillary_number / Ca;
if (force_mag > 1e-3){
Fx *= 1e-3/force_mag; // impose ceiling for stability
Fy *= 1e-3/force_mag;
Fz *= 1e-3/force_mag;
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
color_db->putVector<double>("F",{Fx,Fy,Fz});
}
if ( isSteady ){
MORPH_ADAPT = true;
@@ -846,7 +870,7 @@ void ScaLBL_ColorModel::Run(){
Fx *= capillary_number / Ca;
Fy *= capillary_number / Ca;
Fz *= capillary_number / Ca;
RESCALE_FORCE_COUNT = 1;
if (force_mag > 1e-3){
Fx *= 1e-3/force_mag; // impose ceiling for stability
Fy *= 1e-3/force_mag;
@@ -1451,4 +1475,68 @@ void ScaLBL_ColorModel::WriteDebug(){
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE);
fclose(OUTFILE);
ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
FILE *AFILE;
sprintf(LocalRankFilename,"A.%05i.raw",rank);
AFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,AFILE);
fclose(AFILE);
ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField);
FILE *BFILE;
sprintf(LocalRankFilename,"B.%05i.raw",rank);
BFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,BFILE);
fclose(BFILE);
ScaLBL_Comm->RegularLayout(Map,Pressure,PhaseField);
FILE *PFILE;
sprintf(LocalRankFilename,"Pressure.%05i.raw",rank);
PFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,PFILE);
fclose(PFILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
FILE *VELX_FILE;
sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);
VELX_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,VELX_FILE);
fclose(VELX_FILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
FILE *VELY_FILE;
sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank);
VELY_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,VELY_FILE);
fclose(VELY_FILE);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
FILE *VELZ_FILE;
sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank);
VELZ_FILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,VELZ_FILE);
fclose(VELZ_FILE);
// ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField);
// FILE *CGX_FILE;
// sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank);
// CGX_FILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,CGX_FILE);
// fclose(CGX_FILE);
//
// ScaLBL_Comm->RegularLayout(Map,&ColorGrad[Np],PhaseField);
// FILE *CGY_FILE;
// sprintf(LocalRankFilename,"Gradient_Y.%05i.raw",rank);
// CGY_FILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,CGY_FILE);
// fclose(CGY_FILE);
//
// ScaLBL_Comm->RegularLayout(Map,&ColorGrad[2*Np],PhaseField);
// FILE *CGZ_FILE;
// sprintf(LocalRankFilename,"Gradient_Z.%05i.raw",rank);
// CGZ_FILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,CGZ_FILE);
// fclose(CGZ_FILE);
}

2
models/MRTModel.cpp
View File

@@ -32,7 +32,7 @@ void ScaLBL_MRTModel::ReadParams(string filename){
timestepMax = mrt_db->getScalar<int>( "timestepMax" );
}
if (mrt_db->keyExists( "tolerance" )){
tolerance = mrt_db->getScalar<int>( "timestepMax" );
tolerance = mrt_db->getScalar<double>( "tolerance" );
}
if (mrt_db->keyExists( "tau" )){
tau = mrt_db->getScalar<double>( "tau" );

1
tests/CMakeLists.txt
View File

@@ -69,6 +69,7 @@ ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 )
ADD_LBPM_TEST_1_2_4( testCommunication )
ADD_LBPM_TEST( TestWriter )
ADD_LBPM_TEST( TestDatabase )
ADD_LBPM_PROVISIONAL_TEST( TestMicroCTReader )
IF ( USE_NETCDF )
ADD_LBPM_TEST_PARALLEL( TestNetcdf 8 )
ADD_LBPM_EXECUTABLE( lbpm_uCT_pp )

67
tests/TestMicroCTReader.cpp Normal file
View File

@@ -0,0 +1,67 @@
// Test reading high-resolution files from the microct database
#include "common/MPI_Helpers.h"
#include "common/UnitTest.h"
#include "common/Database.h"
#include "common/Domain.h"
#include "common/ReadMicroCT.h"
#include "IO/Writer.h"
#include <memory>
void testReadMicroCT( const std::string& filename, UnitTest& ut )
{
// Get the domain info
auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" );
// Test reading microCT files
auto data = readMicroCT( *domain_db, MPI_COMM_WORLD );
// Check if we loaded the data correctly
if ( data.size() == domain_db->getVector<size_t>( "n" ) )
ut.passes( "Read data" );
else
ut.passes( "Data size does not match" );
// Write the results to silo
auto n = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
int N[3] = { n[0]*nproc[0], n[1]*nproc[1], n[2]*nproc[2] };
int rank = comm_rank(MPI_COMM_WORLD);
RankInfoStruct rankInfo( rank, nproc[0], nproc[1], nproc[2] );
std::vector<IO::MeshDataStruct> meshData( 1 );
auto Var = std::make_shared<IO::Variable>();
Var->name = "Variable";
Var->type = IO::VariableType::VolumeVariable;
Var->dim = 1;
Var->data.copy( data );
meshData[0].meshName = "grid";
meshData[0].mesh = std::make_shared<IO::DomainMesh>(rankInfo,n[0],n[1],n[2],N[0],N[1],N[2]);
meshData[0].vars.push_back(Var);
IO::writeData( 0, meshData, MPI_COMM_WORLD );
}
int main(int argc, char **argv)
{
// Initialize MPI
MPI_Init(&argc,&argv);
UnitTest ut;
// Run the tests
ASSERT( argc == 2 );
testReadMicroCT( argv[1], ut );
// Print the results
ut.report();
int N_errors = ut.NumFailGlobal();
// Close MPI
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
return N_errors;
}

5
tests/TestPoiseuille.cpp
View File

@@ -34,15 +34,14 @@ void ParallelPlates(ScaLBL_MRTModel &MRT){
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n=k*Nx*Ny+j*Nx+i;
// Initialize distance to +/- 1
MRT.Distance(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
}
}
}
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
if (MRT.rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
CalcDist(MRT.Distance,id_solid,*MRT.Dm);
if (rank == 0) cout << "Domain set." << endl;
if (MRT.rank == 0) cout << "Domain set." << endl;
}
//***************************************************************************************