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

Fixing merge conflicts

Browse Source
This commit is contained in:
James E McClure 2015-11-10 19:36:50 -05:00
commit 596f33318a
11 changed files with 495 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
common/Array.hpp
View File

@ -141,7 +141,7 @@ void Array<TYPE>::resize( const std::vector<size_t>& N )
// Store the old data
const size_t ndim_max = sizeof(d_N)/sizeof(size_t);
std::vector<size_t> N1(ndim_max,1), N2(ndim_max,1);
for (size_t d=0; d<d_ndim; d++)
for (int d=0; d<d_ndim; d++)
N1[d] = d_N[d];
for (size_t d=0; d<N.size(); d++)
N2[d] = N[d];

44
common/Communication.hpp
View File

@ -207,27 +207,27 @@ template<class TYPE1, class TYPE2>
void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
{
PROFILE_START("fillHalo::copy",1);
ASSERT(src.size(0)==nx||src.size(0)==nx+2*ngx);
ASSERT(dst.size(0)==nx||dst.size(0)==nx+2*ngx);
bool src_halo = src.size(0)==nx+2*ngx;
bool dst_halo = dst.size(0)==nx+2*ngx;
ASSERT( (int)src.size(0)==nx || (int)src.size(0)==nx+2*ngx );
ASSERT( (int)dst.size(0)==nx || (int)dst.size(0)==nx+2*ngx );
bool src_halo = (int)src.size(0)==nx+2*ngx;
bool dst_halo = (int)dst.size(0)==nx+2*ngx;
if ( src_halo ) {
ASSERT(src.size(0)==nx+2*ngx);
ASSERT(src.size(1)==ny+2*ngy);
ASSERT(src.size(2)==nz+2*ngz);
ASSERT((int)src.size(0)==nx+2*ngx);
ASSERT((int)src.size(1)==ny+2*ngy);
ASSERT((int)src.size(2)==nz+2*ngz);
} else {
ASSERT(src.size(0)==nx);
ASSERT(src.size(1)==ny);
ASSERT(src.size(2)==nz);
ASSERT((int)src.size(0)==nx);
ASSERT((int)src.size(1)==ny);
ASSERT((int)src.size(2)==nz);
}
if ( dst_halo ) {
ASSERT(dst.size(0)==nx+2*ngx);
ASSERT(dst.size(1)==ny+2*ngy);
ASSERT(dst.size(2)==nz+2*ngz);
ASSERT((int)dst.size(0)==nx+2*ngx);
ASSERT((int)dst.size(1)==ny+2*ngy);
ASSERT((int)dst.size(2)==nz+2*ngz);
} else {
ASSERT(dst.size(0)==nx);
ASSERT(dst.size(1)==ny);
ASSERT(dst.size(2)==nz);
ASSERT((int)dst.size(0)==nx);
ASSERT((int)dst.size(1)==ny);
ASSERT((int)dst.size(2)==nz);
}
if ( src_halo == dst_halo ) {
// Src and dst halos match
@ -235,18 +235,18 @@ void fillHalo<TYPE>::copy( const Array<TYPE1>& src, Array<TYPE2>& dst )
dst(i) = src(i);
} else if ( src_halo && !dst_halo ) {
// Src has halos
for (size_t k=0; k<nz; k++) {
for (size_t j=0; j<ny; j++) {
for (size_t i=0; i<nx; i++) {
for (int k=0; k<nz; k++) {
for (int j=0; j<ny; j++) {
for (int i=0; i<nx; i++) {
dst(i,j,k) = src(i+ngx,j+ngy,k+ngz);
}
}
}
} else if ( !src_halo && dst_halo ) {
// Dst has halos
for (size_t k=0; k<nz; k++) {
for (size_t j=0; j<ny; j++) {
for (size_t i=0; i<nx; i++) {
for (int k=0; k<nz; k++) {
for (int j=0; j<ny; j++) {
for (int i=0; i<nx; i++) {
dst(i+ngx,j+ngy,k+ngz) = src(i,j,k);
}
}

4
common/MPI_Helpers.hpp
View File

@ -27,7 +27,7 @@ void pack( const std::vector<TYPE>& rhs, char *buffer )
size_t size = rhs.size();
memcpy(buffer,&size,sizeof(size_t));
size_t pos = sizeof(size_t);
for (int i=0; i<rhs.size(); i++) {
for (size_t i=0; i<rhs.size(); i++) {
pack(rhs[i],&buffer[pos]);
pos += packsize(rhs[i]);
}
@ -40,7 +40,7 @@ void unpack( std::vector<TYPE>& data, const char *buffer )
data.clear();
data.resize(size);
size_t pos = sizeof(size_t);
for (int i=0; i<data.size(); i++) {
for (size_t i=0; i<data.size(); i++) {
unpack(data[i],&buffer[pos]);
pos += packsize(data[i]);
}

45
common/TwoPhase.cpp
View File

@ -67,7 +67,7 @@ TwoPhase::TwoPhase(Domain &dm):
Jwn(0), Jwn_global(0), Kwn(0), Kwn_global(0), KNwns(0), KNwns_global(0),
KGwns(0), KGwns_global(0), trawn(0), trawn_global(0), trJwn(0), trJwn_global(0),
trRwn(0), trRwn_global(0), nwp_volume_global(0), wp_volume_global(0),
As_global(0), dEs(0), dAwn(0), dAns(0)
As_global(0), dEs(0), dAwn(0), dAns(0), wwndnw(0), wwndnw_global(0)
{
Nx=dm.Nx; Ny=dm.Ny; Nz=dm.Nz;
Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm.nprocx*Dm.nprocy*Dm.nprocz*1.0;
@ -147,8 +147,7 @@ TwoPhase::TwoPhase(Domain &dm):
fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors
fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
fprintf(TIMELOG,"trawn trJwn trRwn Euler Kn Jn An\n"); // trimmed curvature & minkowski measures
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"trawn trJwn trRwn wwndnw Euler Kn Jn An\n"); // trimmed curvature & minkowski measures
}
NWPLOG = fopen("components.NWP.tcat","a+");
@ -175,38 +174,41 @@ TwoPhase::~TwoPhase()
void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, DoubleArray &DistData)
{
/* double factor,temp,value;
double factor,temp,value;
factor=0.5/Beta;
for (int n=0; n<Nx*Ny*Nz; n++){
value = ColorData[n];
if (value > 0.999 ) DistData[n] = 4.0;
else if (value < -0.999 ) DistData[n] = -4.0;
else DistData[n] = factor*log((1.0+value)/(1.0-value));
if (DistData[n] > 1.0) DistData[n] = 1.0;
if (DistData[n] < -1.0) DistData[n] = -1.0;
}
// Initialize to -1,1 (segmentation)
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
value = ColorData(i,j,k);
temp = factor*log((1.0+value)/(1.0-value));
if (temp > 1.0) DistData(i,j,k) = 1.0;
else if (temp < -1.0) DistData(i,j,k) = -1.0;
if (value > 0.8) DistData(i,j,k) = 2.94*factor;
else if (value < -0.8) DistData(i,j,k) = -2.94*factor;
else DistData(i,j,k) = temp;
// Basic threshold
//if (value > 0) DistData(i,j,k) = 1.0;
//else DistData(i,j,k) = -1.0;
}
}
}
SSO(DistData,Dm.id,Dm,10);
*/
SSO(DistData,Dm.id,Dm,40);
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
DistData(i,j,k) += 1.0;
}
}
}
/* for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
DistData(i,j,k) = ColorData(i,j,k);
}
}
}
*/
// for (int n=0; n<Nx*Ny*Nz; n++) DistData[n] = ColorData[n];
}
@ -255,6 +257,7 @@ void TwoPhase::Initialize()
Jwn = Kwn = efawns = 0.0;
trJwn = trawn = trRwn = 0.0;
euler = Jn = An = Kn = 0.0;
wwndnw = 0.0;
}
@ -303,6 +306,8 @@ void TwoPhase::UpdateMeshValues()
{
int i,j,k,n;
//...........................................................................
Dm.CommunicateMeshHalo(SDn);
//...........................................................................
// Compute the gradients of the phase indicator and signed distance fields
pmmc_MeshGradient(SDn,SDn_x,SDn_y,SDn_z,Nx,Ny,Nz);
//...........................................................................
@ -437,6 +442,9 @@ void TwoPhase::ComputeLocal()
// Compute the normal speed of the interface
pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris,
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
for (int p=0; p <n_nw_tris; p++) wwndnw += InterfaceSpeed(p);
}
// wns common curve averages
if (n_local_nws_pts > 0){
@ -446,6 +454,7 @@ void TwoPhase::ComputeLocal()
pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z,
local_nws_pts,i,j,k,n_local_nws_pts);
pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
SDs_z, KNwns_values, KGwns_values, KNwns, KGwns,
nws_pts, n_nws_pts, i, j, k);
@ -1070,6 +1079,7 @@ void TwoPhase::Reduce()
MPI_Allreduce(&KGwns,&KGwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&KNwns,&KNwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&wwndnw,&wwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
// Phase averages
MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
@ -1115,6 +1125,7 @@ void TwoPhase::Reduce()
if (awn_global > 0.0){
Jwn_global /= awn_global;
Kwn_global /= awn_global;
wwndnw_global /= awn_global;
for (i=0; i<3; i++) vawn_global(i) /= awn_global;
for (i=0; i<6; i++) Gwn_global(i) /= awn_global;
}
@ -1174,7 +1185,7 @@ void TwoPhase::PrintAll(int timestep)
Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5)); // orientation of ns interface
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface
fprintf(TIMELOG,"%.5g %.5g %.5g ",trawn_global, trJwn_global, trRwn_global); // Trimmed curvature
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",trawn_global, trJwn_global, trRwn_global, wwndnw_global); // Trimmed curvature
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global); // minkowski measures
fflush(TIMELOG);
}

1
common/TwoPhase.h
View File

@ -95,6 +95,7 @@ public:
double nwp_volume_global; // volume for the non-wetting phase
double wp_volume_global; // volume for the wetting phase
double As_global;
double wwndnw, wwndnw_global;
double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0)
DoubleArray van;
DoubleArray vaw;

6
common/Utilities.cpp
View File

@ -236,9 +236,9 @@ size_t Utilities::getMemoryUsage()
size_t N_bytes = 0;
#if defined(USE_LINUX)
struct mallinfo meminfo = mallinfo();
size_t size_hblkhd = static_cast<size_t>( meminfo.hblkhd );
size_t size_uordblks = static_cast<size_t>( meminfo.uordblks );
N_bytes = static_cast<size_t>( size_hblkhd + size_uordblks );
size_t size_hblkhd = static_cast<unsigned int>( meminfo.hblkhd );
size_t size_uordblks = static_cast<unsigned int>( meminfo.uordblks );
N_bytes = size_hblkhd + size_uordblks;
#elif defined(USE_MAC)
struct task_basic_info t_info;
mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT;

20
common/pmmc.h
View File

@ -4317,12 +4317,12 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
//--------------------------------------------------------------------------------------------------------
inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray &P_y, DoubleArray &P_z,
DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
DoubleArray &SurfaceVector, DoubleArray &SurfaceValues, DoubleArray &AvgVel,
DoubleArray &SurfaceVector, DoubleArray &AvgSpeed, DoubleArray &AvgVel,
int i, int j, int k, int npts, int ntris)
{
Point A,B,C,P;
double x,y,z;
double s,s1,s2,s3,temp;
double s,s1,s2,s3,area;
double norm, zeta;
TriLinPoly Px,Py,Pz,Pt;
@ -4342,23 +4342,27 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
temp = s*(s-s1)*(s-s2)*(s-s3);
area = sqrt(s*(s-s1)*(s-s2)*(s-s3));
// Compute the centroid P
P.x = 0.33333333333333333*(A.x+B.x+C.x);
P.y = 0.33333333333333333*(A.y+B.y+C.y);
P.z = 0.33333333333333333*(A.z+B.z+C.z);
if (temp > 0.0){
if (area > 0.0){
x = Px.eval(P);
y = Py.eval(P);
z = Pz.eval(P);
norm = sqrt(x*x+y*y+z*z);
if (norm==0.0) norm=1.0;
// Compute the interface speed from time derivative and gradient (Level Set Equation)
zeta = -Pt.eval(P) / norm;
temp = sqrt(temp)/norm;
AvgVel(0) += temp*zeta*x;
AvgVel(1) += temp*zeta*y;
AvgVel(2) += temp*zeta*z;
//temp = sqrt(temp)/norm; <--- what was I thinking with this? (James)
// Compute the average
AvgVel(0) += area*zeta*x;
AvgVel(1) += area*zeta*y;
AvgVel(2) += area*zeta*z;
AvgSpeed(r) = zeta*area;
}
}
//.............................................................................

1
tests/CMakeLists.txt
View File

@ -12,6 +12,7 @@ ADD_LBPM_EXECUTABLE( lbpm_BlobAnalysis )
ADD_LBPM_EXECUTABLE( TestBubble )
ADD_LBPM_EXECUTABLE( BasicSimulator )
ADD_LBPM_EXECUTABLE( ComponentLabel )
ADD_LBPM_EXECUTABLE( ColorToBinary )
ADD_LBPM_EXECUTABLE( BlobAnalysis )
ADD_LBPM_EXECUTABLE( BlobIdentify )
ADD_LBPM_EXECUTABLE( BlobIdentifyParallel )

261
tests/ColorToBinary.cpp Normal file
View File

@ -0,0 +1,261 @@
// Sequential component labeling for two phase systems
// Reads parallel simulation data and performs connectivity analysis
// and averaging on a blob-by-blob basis
// James E. McClure 2015
#include <iostream>
#include <math.h>
#include "analysis/analysis.h"
#include "common/TwoPhase.h"
using namespace std;
inline void ReadFromRank(char *FILENAME, DoubleArray &Phase, int nx, int ny, int nz, int iproc, int jproc, int kproc)
{
int i,j,k,q,n,N;
int iglobal,jglobal,kglobal;
double value;
double denA,denB;
double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
double f10,f11,f12,f13,f14,f15,f16,f17,f18;
double vx,vy,vz;
N = nx*ny*nz;
double *Den, *DistEven, *DistOdd;
Den = new double[2*N];
DistEven = new double[10*N];
DistOdd = new double[9*N];
ifstream File(FILENAME,ios::binary);
for (n=0; n<N; n++){
// Write the two density values
File.read((char*) &value, sizeof(value));
Den[n] = value;
// if (n== 66276) printf("Density a = %f \n",value);
File.read((char*) &value, sizeof(value));
Den[N+n] = value;
// if (n== 66276) printf("Density b = %f \n",value);
// Read the even distributions
for (q=0; q<10; q++){
File.read((char*) &value, sizeof(value));
DistEven[q*N+n] = value;
}
// Read the odd distributions
for (q=0; q<9; q++){
File.read((char*) &value, sizeof(value));
DistOdd[q*N+n] = value;
}
}
File.close();
// Compute the phase field, pressure and velocity
for (k=1; k<nz-1; k++){
for (j=1; j<ny-1; j++){
for (i=1; i<nz-1; i++){
//........................................................................
n = k*nx*ny+j*nx+i;
//........................................................................
denA = Den[n];
denB = Den[N+n];
//........................................................................
f0 = DistEven[n];
f2 = DistEven[N+n];
f4 = DistEven[2*N+n];
f6 = DistEven[3*N+n];
f8 = DistEven[4*N+n];
f10 = DistEven[5*N+n];
f12 = DistEven[6*N+n];
f14 = DistEven[7*N+n];
f16 = DistEven[8*N+n];
f18 = DistEven[9*N+n];
//........................................................................
f1 = DistOdd[n];
f3 = DistOdd[1*N+n];
f5 = DistOdd[2*N+n];
f7 = DistOdd[3*N+n];
f9 = DistOdd[4*N+n];
f11 = DistOdd[5*N+n];
f13 = DistOdd[6*N+n];
f15 = DistOdd[7*N+n];
f17 = DistOdd[8*N+n];
//........................................................................
//.................Compute the pressure....................................
value = 0.3333333333333333*(f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17);
//........................................................................
//.................Compute the velocity...................................
vx = f1-f2+f7-f8+f9-f10+f11-f12+f13-f14;
vy = f3-f4+f7-f8-f9+f10+f15-f16+f17-f18;
vz = f5-f6+f11-f12-f13+f14+f15-f16-f17+f18;
//........................................................................
// save values in global arrays
//........................................................................
iglobal = iproc*(nx-2)+i;
jglobal = jproc*(ny-2)+j;
kglobal = kproc*(nz-2)+k;
//........................................................................
Phase(iglobal,jglobal,kglobal) = (denA-denB)/(denA+denB);
//Pressure(iglobal,jglobal,kglobal) = value;
//Vel_x(iglobal,jglobal,kglobal) = vx;
//Vel_y(iglobal,jglobal,kglobal) = vy;
//Vel_z(iglobal,jglobal,kglobal) = vz;
//........................................................................
}
}
}
delete Den;
delete DistEven;
delete DistOdd;
}
int main(int argc, char **argv)
{
// Initialize MPI
int rank,nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
printf("----------------------------------------------------------\n");
printf("Creating single Binary file from restart (8-bit integer)\n");
printf("ID=0 (solid) ID=1 (non-wetting), ID=2 (wetting) \n");
printf("----------------------------------------------------------\n");
if (nprocs != 1) INSIST(nprocs == 1,"Error: ComponentLabel --serial case!");
//.......................................................................
int nprocx,nprocy,nprocz;
int Nx, Ny, Nz;
int nx,ny,nz;
int nspheres;
double Lx,Ly,Lz;
//.......................................................................
int i,j,k,n;
int iproc,jproc,kproc;
//.......................................................................
// Reading the domain information file
//.......................................................................
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> nx;
domain >> ny;
domain >> nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
//.......................................................................
nx+=2;
ny+=2;
nz+=2;
nprocs = nprocx*nprocy*nprocz;
printf("Number of MPI ranks: %i \n", nprocs);
int BoundaryCondition=0;
Nx = (nx-2)*nprocx;
Ny = (ny-2)*nprocy;
Nz = (nz-2)*nprocz;
Domain Dm(Nx,Ny,Nz,rank,1,1,1,Lx,Ly,Lz,BoundaryCondition);
Nx+=2; Ny+=2; Nz+=2;
printf("Full domain size: %i x %i x %i \n", Nx,Ny,Nz);
DoubleArray Phase(Nx,Ny,Nz);
DoubleArray SignDist(Nx,Ny,Nz);
// Filenames used
char LocalRankString[8];
char LocalRankFilename[40];
char BaseFilename[20];
int proc,iglobal,kglobal,jglobal;
double * Temp;
Temp = new double[nx*ny*nz];
// read the files and populate main arrays
for ( kproc=0; kproc<nprocz; kproc++){
for ( jproc=0; jproc<nprocy; jproc++){
for ( iproc=0; iproc<nprocx; iproc++){
proc = kproc*nprocx*nprocy + jproc*nprocx + iproc;
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);
for (k=1; k<nz-1; k++){
for (j=1; j<ny-1; j++){
for (i=1; i<nz-1; i++){
//........................................................................
n = k*nx*ny+j*nx+i;
//........................................................................
iglobal = iproc*(nx-2)+i;
jglobal = jproc*(ny-2)+j;
kglobal = kproc*(nz-2)+k;
//........................................................................
SignDist(iglobal,jglobal,kglobal) = Temp[n];
//........................................................................
}
}
}
sprintf(LocalRankFilename,"%s%s","Restart.",LocalRankString);
ReadFromRank(LocalRankFilename,Phase,nx,ny,nz,iproc,jproc,kproc);
}
}
}
printf("Read %i ranks of %s \n",nprocs,BaseFilename);
delete Temp;
// Initializing the blob ID
char *PhaseID;
PhaseID = new char (Nx*Ny*Nz);
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
n = k*Nx*Ny+j*Nx+i;
if (SignDist(i,j,k) < 0.0){
// Solid phase
PhaseID[n] = 0;
}
else if (Phase(i,j,k) < 0.0){
// wetting phase
PhaseID[n] = 2;
}
else {
// non-wetting phase
PhaseID[n] = 1;
}
}
}
}
FILE *OUTFILE;
OUTFILE = fopen("ID.dat","wb");
fwrite(Dm.id,1,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
OUTFILE = fopen("Phase.dat","wb");
fwrite(Phase.get(),8,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
OUTFILE = fopen("SignDist.dat","wb");
fwrite(SignDist.get(),8,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
// ****************************************************
MPI_Barrier(comm);
MPI_Finalize();
// ****************************************************
}

115
tests/lbpm_color_simulator.cpp
View File

@ -290,6 +290,7 @@ int main(int argc, char **argv)
if (BoundaryCondition==0) printf("Periodic boundary conditions will applied \n");
if (BoundaryCondition==1) printf("Pressure boundary conditions will be applied \n");
if (BoundaryCondition==2) printf("Velocity boundary conditions will be applied \n");
if (BoundaryCondition==3) printf("Dynamic pressure boundary conditions will be applied \n");
if (InitialCondition==0) printf("Initial conditions assigned from phase ID file \n");
if (InitialCondition==1) printf("Initial conditions assigned from restart file \n");
printf("********************************************************\n");
@ -297,10 +298,12 @@ int main(int argc, char **argv)
// Initialized domain and averaging framework for Two-Phase Flow
bool pBC,velBC;
if (BoundaryCondition==1) pBC=true;
if (BoundaryCondition==1 || BoundaryCondition==3)
pBC=true;
else pBC=false;
if (BoundaryCondition==2) velBC=true;
else velBC=false;
bool Restart;
if (InitialCondition==1) Restart=true;
else Restart=false;
@ -575,7 +578,7 @@ int main(int argc, char **argv)
MPI_Barrier(comm);
//.......................................................................
// Once phase has been initialized, map solid to account for 'smeared' interface
for (i=0; i<N; i++) Averages->SDs(i) -= (1.0); //
// for (i=0; i<N; i++) Averages->SDs(i) -= (1.0); //
//.......................................................................
// Finalize setup for averaging domain
//Averages->SetupCubes(Dm);
@ -624,6 +627,26 @@ int main(int argc, char **argv)
SetPhiSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
}
// Set dynamic pressure boundary conditions
double dp, slope;
if (BoundaryCondition==3){
slope = (dout-din)/timestepMax;
dp = din;
if (rank==0) printf("Change in pressure / time =%f \n",slope);
// set the initial value
din = 1.0+0.5*dp;
dout = 1.0-0.5*dp;
// set the initial boundary conditions
if (Dm.kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (Dm.kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
}
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
ComputeVelocityD3Q19(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
@ -684,6 +707,37 @@ int main(int argc, char **argv)
ThreadPool::setProcessAffinity(procs);
}
ThreadPool tpool(N_threads);
// Create the MeshDataStruct
fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
std::vector<IO::MeshDataStruct> meshData(1);
meshData[0].meshName = "domain";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
std::shared_ptr<IO::Variable> PressVar( new IO::Variable() );
std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
PhaseVar->name = "phase";
PhaseVar->type = IO::VolumeVariable;
PhaseVar->dim = 1;
PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(PhaseVar);
PressVar->name = "Pressure";
PressVar->type = IO::VolumeVariable;
PressVar->dim = 1;
PressVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(PressVar);
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VolumeVariable;
SignDistVar->dim = 1;
SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(SignDistVar);
BlobIDVar->name = "BlobID";
BlobIDVar->type = IO::VolumeVariable;
BlobIDVar->dim = 1;
BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(BlobIDVar);
//************ MAIN ITERATION LOOP ***************************************/
PROFILE_START("Loop");
int timestep = -1;
@ -692,6 +746,7 @@ int main(int argc, char **argv)
writeIDMap(ID_map_struct(),0,id_map_filename);
AnalysisWaitIdStruct work_ids;
while (timestep < timestepMax && err > tol ) {
if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
PROFILE_START("Update");
//*************************************************************************
@ -787,6 +842,23 @@ int main(int argc, char **argv)
//ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
SetPhiSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
}
if (BoundaryCondition==3){
// Increase the pressure difference
dp += slope;
din = 1.0+0.5*dp;
dout = 1.0-0.5*dp;
// set the initial boundary conditions
if (Dm.kproc == 0) {
PressureBC_inlet(f_even,f_odd,din,Nx,Ny,Nz);
ColorBC_inlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
if (Dm.kproc == nprocz-1){
PressureBC_outlet(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
ColorBC_outlet(Phi,Den,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
}
}
//...................................................................................
MPI_Barrier(comm);
@ -798,8 +870,11 @@ int main(int argc, char **argv)
// Run the analysis, blob identification, and write restart files
run_analysis(timestep,RESTART_INTERVAL,rank_info,*Averages,last_ids,last_index,last_id_map,
Nx,Ny,Nz,pBC,beta,err,Phi,Pressure,Velocity,ID,f_even,f_odd,Den,
LocalRestartFile,tpool,work_ids);
LocalRestartFile,meshData,fillData,tpool,work_ids);
// Save the timers
if ( timestep%50==0 )
PROFILE_SAVE("lbpm_color_simulator",1);
}
tpool.wait_pool_finished();
PROFILE_STOP("Loop");
@ -833,40 +908,6 @@ int main(int argc, char **argv)
DeviceBarrier();
CopyToHost(Averages->Phase.get(),Phi,N*sizeof(double));
*/
// Create the MeshDataStruct
fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
std::vector<IO::MeshDataStruct> meshData(1);
meshData[0].meshName = "domain";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
std::shared_ptr<IO::Variable> PressVar( new IO::Variable() );
std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
PhaseVar->name = "phase";
PhaseVar->type = IO::VolumeVariable;
PhaseVar->dim = 1;
PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(PhaseVar);
PressVar->name = "Pressure";
PressVar->type = IO::VolumeVariable;
PressVar->dim = 1;
PressVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(PressVar);
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VolumeVariable;
SignDistVar->dim = 1;
SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(SignDistVar);
BlobIDVar->name = "BlobID";
BlobIDVar->type = IO::VolumeVariable;
BlobIDVar->dim = 1;
BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
meshData[0].vars.push_back(BlobIDVar);
fillData.copy(Averages->SDn,PhaseVar->data);
fillData.copy(Averages->SDs,SignDistVar->data);
fillData.copy(Averages->Label_NWP,BlobIDVar->data);
IO::writeData( 0, meshData, 2, comm );
/* Averages->WriteSurfaces(0);

91
tests/lbpm_color_simulator.h
View File

@ -1,17 +1,22 @@
// Run the analysis, blob identification, and write restart files
#include "common/Array.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "IO/MeshDatabase.h"
//#define ANALYSIS_INTERVAL 6
#define ANALYSIS_INTERVAL 1000
#define BLOBID_INTERVAL 1000
#define BLOBID_INTERVAL 250
enum AnalysisType{ AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02,
CopyAverages=0x04, CalcDist=0x08, CreateRestart=0x10 };
CopyAverages=0x04, CalcDist=0x08, CreateRestart=0x10, WriteVis=0x20 };
// Structure used to store ids
struct AnalysisWaitIdStruct {
ThreadPool::thread_id_t blobID;
ThreadPool::thread_id_t analysis;
ThreadPool::thread_id_t vis;
ThreadPool::thread_id_t restart;
};
@ -28,7 +33,6 @@ public:
PROFILE_START("Save Checkpoint",1);
WriteCheckpoint(filename,cDen.get(),cDistEven.get(),cDistOdd.get(),N);
PROFILE_STOP("Save Checkpoint",1);
PROFILE_SAVE("lbpm_color_simulator",1);
ThreadPool::WorkItem::d_state = 2; // Change state to finished
};
private:
@ -123,6 +127,44 @@ private:
};
// Helper class to write the vis file from a thread
class WriteVisWorkItem: public ThreadPool::WorkItem
{
public:
WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
TwoPhase& Avgerages_, fillHalo<double>& fillData_ ):
timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_) {}
virtual void run() {
ThreadPool::WorkItem::d_state = 1; // Change state to in progress
PROFILE_START("Save Vis",1);
ASSERT(visData[0].vars[0]->name=="phase");
ASSERT(visData[0].vars[1]->name=="Pressure");
ASSERT(visData[0].vars[2]->name=="SignDist");
ASSERT(visData[0].vars[3]->name=="BlobID");
Array<double>& PhaseData = visData[0].vars[0]->data;
Array<double>& PressData = visData[0].vars[1]->data;
Array<double>& SignData = visData[0].vars[2]->data;
Array<double>& BlobData = visData[0].vars[3]->data;
fillData.copy(Averages.SDn,PhaseData);
fillData.copy(Averages.Press,PressData);
fillData.copy(Averages.SDs,SignData);
fillData.copy(Averages.Label_NWP,BlobData);
MPI_Comm newcomm;
MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
IO::writeData( timestep, visData, 2, newcomm );
MPI_Comm_free(&newcomm);
PROFILE_STOP("Save Vis",1);
ThreadPool::WorkItem::d_state = 2; // Change state to finished
};
private:
WriteVisWorkItem();
int timestep;
std::vector<IO::MeshDataStruct>& visData;
TwoPhase& Averages;
fillHalo<double>& fillData;
};
// Helper class to run the analysis from within a thread
// Note: Averages will be modified after the constructor is called
class AnalysisWorkItem: public ThreadPool::WorkItem
@ -170,6 +212,8 @@ private:
double beta;
};
// Function to start the analysis
void run_analysis( int timestep, int restart_interval,
const RankInfoStruct& rank_info, TwoPhase& Averages,
@ -177,7 +221,8 @@ void run_analysis( int timestep, int restart_interval,
int Nx, int Ny, int Nz, bool pBC, double beta, double err,
const double *Phi, double *Pressure, const double *Velocity,
const char *ID, const double *f_even, const double *f_odd, const double *Den,
const char *LocalRestartFile, ThreadPool& tpool, AnalysisWaitIdStruct& wait )
const char *LocalRestartFile, std::vector<IO::MeshDataStruct>& visData, fillHalo<double>& fillData,
ThreadPool& tpool, AnalysisWaitIdStruct& wait )
{
int N = Nx*Ny*Nz;
@ -191,7 +236,7 @@ void run_analysis( int timestep, int restart_interval,
// Identify blobs and update global ids in time
type = static_cast<AnalysisType>( type | IdentifyBlobs );
}
/* #ifdef USE_CUDA
#ifdef USE_CUDA
if ( tpool.getQueueSize()<=3 && tpool.getNumThreads()>0 && timestep%50==0 ) {
// Keep a few blob identifications queued up to keep the processors busy,
// allowing us to track the blobs as fast as possible
@ -199,7 +244,7 @@ void run_analysis( int timestep, int restart_interval,
type = static_cast<AnalysisType>( type | IdentifyBlobs );
}
#endif
*/
if ( timestep%ANALYSIS_INTERVAL == 0 ) {
// Copy the averages to the CPU (and identify blobs)
type = static_cast<AnalysisType>( type | CopyAverages );
@ -213,6 +258,12 @@ void run_analysis( int timestep, int restart_interval,
// Write the restart file
type = static_cast<AnalysisType>( type | CreateRestart );
}
if (timestep%restart_interval == 0) {
// Write the visualization data
type = static_cast<AnalysisType>( type | WriteVis );
type = static_cast<AnalysisType>( type | CopyAverages );
type = static_cast<AnalysisType>( type | IdentifyBlobs );
}
// Return if we are not doing anything
if ( type == AnalyzeNone )
@ -238,14 +289,22 @@ void run_analysis( int timestep, int restart_interval,
}
if ( (type&CopyAverages) != 0 ) {
// Copy the members of Averages to the cpu (phase was copied above)
// Wait
PROFILE_START("Copy-Pressure",1);
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
memcpy(Averages.Phase.get(),phase->get(),N*sizeof(double));
DeviceBarrier();
PROFILE_STOP("Copy-Pressure",1);
PROFILE_START("Copy-Wait",1);
tpool.wait(wait.analysis);
tpool.wait(wait.vis); // Make sure we are done using analysis before modifying
PROFILE_STOP("Copy-Wait",1);
PROFILE_START("Copy-Averages",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));
PROFILE_STOP("Copy-Averages",1);
}
std::shared_ptr<double> cDen, cDistEven, cDistOdd;
if ( (type&CreateRestart) != 0 ) {
@ -282,6 +341,7 @@ void run_analysis( int timestep, int restart_interval,
type,timestep,Averages,last_index,last_id_map,beta);
work->add_dependency(wait.blobID);
work->add_dependency(wait.analysis);
work->add_dependency(wait.vis); // Make sure we are done using analysis before modifying
wait.analysis = tpool.add_work(work);
}
@ -295,12 +355,27 @@ void run_analysis( int timestep, int restart_interval,
} else {
// Not clear yet
}
// Wait for previous restart files to finish writing (not necessary, but helps to ensure memory usage is limited)
tpool.wait(wait.restart);
// Write the restart file (using a seperate thread)
WriteRestartWorkItem *work = new WriteRestartWorkItem(LocalRestartFile,cDen,cDistEven,cDistOdd,N);
work->add_dependency(wait.restart);
wait.restart = tpool.add_work(work);
}
// Save the results for visualization
if ( (type&CreateRestart) != 0 ) {
// Wait for previous restart files to finish writing (not necessary, but helps to ensure memory usage is limited)
tpool.wait(wait.vis);
// Write the vis files
ThreadPool::WorkItem *work = new WriteVisWorkItem( timestep, visData, Averages, fillData );
work->add_dependency(wait.blobID);
work->add_dependency(wait.analysis);
work->add_dependency(wait.vis);
wait.vis = tpool.add_work(work);
}
PROFILE_STOP("start_analysis");
}