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:
JamesEMcclure
2019-08-07 16:38:12 -04:00
parent 04648451d0 6001b6a440
commit e049fa4482
12 changed files with 906 additions and 111 deletions
Download Patch File Download Diff File Expand all files Collapse all files

159
analysis/SubPhase.cpp
View File

@@ -232,6 +232,19 @@ void SubPhase::Basic(){
gwb.p=sumReduce( Dm->Comm, wb.p) / count_w;
gnb.p=sumReduce( Dm->Comm, nb.p) / count_n;
// check for NaN
bool err=false;
if (gwb.V != gwb.V) err=true;
if (gnb.V != gnb.V) err=true;
if (gwb.p != gwb.p) err=true;
if (gnb.p != gnb.p) err=true;
if (gwb.Px != gwb.Px) err=true;
if (gwb.Py != gwb.Py) err=true;
if (gwb.Pz != gwb.Pz) err=true;
if (gnb.Px != gnb.Px) err=true;
if (gnb.Py != gnb.Py) err=true;
if (gnb.Pz != gnb.Pz) err=true;
if (Dm->rank() == 0){
double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
double dir_x = Fx/force_mag;
@@ -257,7 +270,11 @@ void SubPhase::Basic(){
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*h*h*water_flow_rate,h*h*h*not_water_flow_rate, gwb.p, gnb.p);
fflush(TIMELOG);
}
if (err==true){
// exception if simulation produceds NaN
printf("SubPhase.cpp: NaN encountered, may need to check simulation parameters \n");
}
ASSERT(err==false);
}
inline void InterfaceTransportMeasures( double beta, double rA, double rB, double nA, double nB,
@@ -492,6 +509,10 @@ void SubPhase::Full(){
double vol_wc_bulk = 0.0;
double vol_nd_bulk = 0.0;
double vol_wd_bulk = 0.0;
double count_wc = 0.0;
double count_nc = 0.0;
double count_wd = 0.0;
double count_nd = 0.0;
for (k=kmin; k<kmax; k++){
for (j=jmin; j<Ny-1; j++){
for (i=imin; i<Nx-1; i++){
@@ -515,6 +536,27 @@ void SubPhase::Full(){
InterfaceTransportMeasures( beta, rho_w, rho_n, nA, nB, nx, ny, nz, ux, uy, uz, iwn);
}
else if ( phi > 0.0){
if (morph_n->label(i,j,k) > 0 ){
count_nd += 1.0;
nd.p += Pressure(n);
}
else{
count_nc += 1.0;
nc.p += Pressure(n);
}
}
else{
// water region
if (morph_w->label(i,j,k) > 0 ){
count_wd += 1.0;
wd.p += Pressure(n);
}
else{
count_wc += 1.0;
wc.p += Pressure(n);
}
}
if ( phi > 0.0){
if (morph_n->label(i,j,k) > 0 ){
vol_nd_bulk += 1.0;
nd.M += nA*rho_n;
@@ -522,7 +564,6 @@ void SubPhase::Full(){
nd.Py += nA*rho_n*uy;
nd.Pz += nA*rho_n*uz;
nd.K += nA*rho_n*(ux*ux + uy*uy + uz*uz);
nd.p += Pressure(n);
}
else{
vol_nc_bulk += 1.0;
@@ -531,7 +572,6 @@ void SubPhase::Full(){
nc.Py += nA*rho_n*uy;
nc.Pz += nA*rho_n*uz;
nc.K += nA*rho_n*(ux*ux + uy*uy + uz*uz);
nc.p += Pressure(n);
}
}
else{
@@ -543,7 +583,6 @@ void SubPhase::Full(){
wd.Py += nB*rho_w*uy;
wd.Pz += nB*rho_w*uz;
wd.K += nB*rho_w*(ux*ux + uy*uy + uz*uz);
wd.p += Pressure(n);
}
else{
vol_wc_bulk += 1.0;
@@ -552,40 +591,44 @@ void SubPhase::Full(){
wc.Py += nB*rho_w*uy;
wc.Pz += nB*rho_w*uz;
wc.K += nB*rho_w*(ux*ux + uy*uy + uz*uz);
wc.p += Pressure(n);
}
}
}
}
}
}
count_wc=sumReduce( Dm->Comm, count_wc);
count_nc=sumReduce( Dm->Comm, count_nc);
count_wd=sumReduce( Dm->Comm, count_wd);
count_nd=sumReduce( Dm->Comm, count_nd);
gnd.p=sumReduce( Dm->Comm, nd.p) / count_nd;
gwd.p=sumReduce( Dm->Comm, wd.p) / count_wd;
gnc.p=sumReduce( Dm->Comm, nc.p) / count_nc;
gwc.p=sumReduce( Dm->Comm, wc.p) / count_wc;
gnd.M=sumReduce( Dm->Comm, nd.M);
gnd.Px=sumReduce( Dm->Comm, nd.Px);
gnd.Py=sumReduce( Dm->Comm, nd.Py);
gnd.Pz=sumReduce( Dm->Comm, nd.Pz);
gnd.K=sumReduce( Dm->Comm, nd.K);
gnd.p=sumReduce( Dm->Comm, nd.p);
gwd.M=sumReduce( Dm->Comm, wd.M);
gwd.Px=sumReduce( Dm->Comm, wd.Px);
gwd.Py=sumReduce( Dm->Comm, wd.Py);
gwd.Pz=sumReduce( Dm->Comm, wd.Pz);
gwd.K=sumReduce( Dm->Comm, wd.K);
gwd.p=sumReduce( Dm->Comm, wd.p);
gnc.M=sumReduce( Dm->Comm, nc.M);
gnc.Px=sumReduce( Dm->Comm, nc.Px);
gnc.Py=sumReduce( Dm->Comm, nc.Py);
gnc.Pz=sumReduce( Dm->Comm, nc.Pz);
gnc.K=sumReduce( Dm->Comm, nc.K);
gnc.p=sumReduce( Dm->Comm, nc.p);
gwc.M=sumReduce( Dm->Comm, wc.M);
gwc.Px=sumReduce( Dm->Comm, wc.Px);
gwc.Py=sumReduce( Dm->Comm, wc.Py);
gwc.Pz=sumReduce( Dm->Comm, wc.Pz);
gwc.K=sumReduce( Dm->Comm, wc.K);
gwc.p=sumReduce( Dm->Comm, wc.p);
giwn.Mn=sumReduce( Dm->Comm, iwn.Mn);
giwn.Pnx=sumReduce( Dm->Comm, iwn.Pnx);
@@ -624,4 +667,102 @@ void SubPhase::Full(){
}
void SubPhase::AggregateLabels(char *FILENAME){
int nx = Dm->Nx;
int ny = Dm->Ny;
int nz = Dm->Nz;
int npx = Dm->nprocx();
int npy = Dm->nprocy();
int npz = Dm->nprocz();
int ipx = Dm->iproc();
int ipy = Dm->jproc();
int ipz = Dm->kproc();
int nprocs = Dm->nprocx()*Dm->nprocy()*Dm->nprocz();
int full_nx = npx*(nx-2);
int full_ny = npy*(ny-2);
int full_nz = npz*(nz-2);
int local_size = (nx-2)*(ny-2)*(nz-2);
long int full_size = long(full_nx)*long(full_ny)*long(full_nz);
signed char *LocalID;
LocalID = new signed char [local_size];
//printf("aggregate labels: local size=%i, global size = %i",local_size, full_size);
// assign the ID for the local sub-region
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
int n = k*nx*ny+j*nx+i;
signed char local_id_val = Dm->id[n];
if (local_id_val > 0){
double value = Phi(i,j,k);
if (value > 0.0) local_id_val = 1;
else local_id_val = 2;
}
LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val;
}
}
}
MPI_Barrier(Dm->Comm);
// populate the FullID
if (Dm->rank() == 0){
signed char *FullID;
FullID = new signed char [full_size];
// first handle local ID for rank 0
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
int x = i-1;
int y = j-1;
int z = k-1;
int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
int n_full = z*full_nx*full_ny + y*full_nx + x;
FullID[n_full] = LocalID[n_local];
}
}
}
// next get the local ID from the other ranks
for (int rnk = 1; rnk<nprocs; rnk++){
ipz = rnk / (npx*npy);
ipy = (rnk - ipz*npx*npy) / npx;
ipx = (rnk - ipz*npx*npy - ipy*npx);
//printf("ipx=%i ipy=%i ipz=%i\n", ipx, ipy, ipz);
int tag = 15+rnk;
MPI_Recv(LocalID,local_size,MPI_CHAR,rnk,tag,Dm->Comm,MPI_STATUS_IGNORE);
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
int x = i-1 + ipx*(nx-2);
int y = j-1 + ipy*(ny-2);
int z = k-1 + ipz*(nz-2);
int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
int n_full = z*full_nx*full_ny + y*full_nx + x;
FullID[n_full] = LocalID[n_local];
}
}
}
}
// write the output
FILE *OUTFILE;
OUTFILE = fopen(FILENAME,"wb");
fwrite(FullID,1,full_size,OUTFILE);
fclose(OUTFILE);
}
else{
// send LocalID to rank=0
int tag = 15+ Dm->rank();
int dstrank = 0;
MPI_Send(LocalID,local_size,MPI_CHAR,dstrank,tag,Dm->Comm);
}
MPI_Barrier(Dm->Comm);
}

1
analysis/SubPhase.h
View File

@@ -101,6 +101,7 @@ public:
void Basic();
void Full();
void Write(int time);
void AggregateLabels(char *FILENAME);
private:
FILE *TIMELOG;

2
analysis/dcel.cpp
View File

@@ -312,6 +312,7 @@ double DECL::EdgeAngle(int edge)
}
if (dotprod > 1.f) dotprod=1.f;
if (dotprod < -1.f) dotprod=-1.f;
angle = acos(dotprod);
/* project onto plane of cube face also works
W = U - dotprod*V;
@@ -344,6 +345,7 @@ double DECL::EdgeAngle(int edge)
// concave
angle = -angle;
}
if (angle != angle) angle = 0.0;
//printf("angle=%f,dot=%f (Edge=%i, twin=%i): P={%f, %f, %f}, Q={%f, %f, %f} U={%f, %f, %f}, V={%f, %f, %f}\n",angle,dotprod,edge,halfedge.twin(edge),P.x,P.y,P.z,Q.x,Q.y,Q.z,U.x,U.y,U.z,V.x,V.y,V.z);
return angle;
}

4
analysis/runAnalysis.cpp
View File

@@ -247,6 +247,9 @@ public:
fillData.copy(Averages.Vel_z,VelzData);
fillData.copy(Averages.morph_n->label,BlobData);
IO::writeData( timestep, visData, comm.comm );
char CurrentIDFilename[40];
sprintf(CurrentIDFilename,"id_t%d.raw",timestep);
Averages.AggregateLabels(CurrentIDFilename);
PROFILE_STOP("Save Vis",1);
};
private:
@@ -1020,4 +1023,3 @@ void runAnalysis::WriteVisData( int timestep, SubPhase &Averages, const double *
PROFILE_STOP("write vis");
}

434
common/Domain.cpp
View File

@@ -128,6 +128,59 @@ Domain::Domain( std::shared_ptr<Database> db, MPI_Comm Communicator):
rank_info = RankInfoStruct( myrank, rank_info.nx, rank_info.ny, rank_info.nz );
MPI_Barrier(Comm);
}
Domain::~Domain()
{
// Free sendList
delete [] sendList_x; delete [] sendList_y; delete [] sendList_z;
delete [] sendList_X; delete [] sendList_Y; delete [] sendList_Z;
delete [] sendList_xy; delete [] sendList_yz; delete [] sendList_xz;
delete [] sendList_Xy; delete [] sendList_Yz; delete [] sendList_xZ;
delete [] sendList_xY; delete [] sendList_yZ; delete [] sendList_Xz;
delete [] sendList_XY; delete [] sendList_YZ; delete [] sendList_XZ;
// Free sendBuf
delete [] sendBuf_x; delete [] sendBuf_y; delete [] sendBuf_z;
delete [] sendBuf_X; delete [] sendBuf_Y; delete [] sendBuf_Z;
delete [] sendBuf_xy; delete [] sendBuf_yz; delete [] sendBuf_xz;
delete [] sendBuf_Xy; delete [] sendBuf_Yz; delete [] sendBuf_xZ;
delete [] sendBuf_xY; delete [] sendBuf_yZ; delete [] sendBuf_Xz;
delete [] sendBuf_XY; delete [] sendBuf_YZ; delete [] sendBuf_XZ;
// Free recvList
delete [] recvList_x; delete [] recvList_y; delete [] recvList_z;
delete [] recvList_X; delete [] recvList_Y; delete [] recvList_Z;
delete [] recvList_xy; delete [] recvList_yz; delete [] recvList_xz;
delete [] recvList_Xy; delete [] recvList_Yz; delete [] recvList_xZ;
delete [] recvList_xY; delete [] recvList_yZ; delete [] recvList_Xz;
delete [] recvList_XY; delete [] recvList_YZ; delete [] recvList_XZ;
// Free recvBuf
delete [] recvBuf_x; delete [] recvBuf_y; delete [] recvBuf_z;
delete [] recvBuf_X; delete [] recvBuf_Y; delete [] recvBuf_Z;
delete [] recvBuf_xy; delete [] recvBuf_yz; delete [] recvBuf_xz;
delete [] recvBuf_Xy; delete [] recvBuf_Yz; delete [] recvBuf_xZ;
delete [] recvBuf_xY; delete [] recvBuf_yZ; delete [] recvBuf_Xz;
delete [] recvBuf_XY; delete [] recvBuf_YZ; delete [] recvBuf_XZ;
// Free sendData
delete [] sendData_x; delete [] sendData_y; delete [] sendData_z;
delete [] sendData_X; delete [] sendData_Y; delete [] sendData_Z;
delete [] sendData_xy; delete [] sendData_xY; delete [] sendData_Xy;
delete [] sendData_XY; delete [] sendData_xz; delete [] sendData_xZ;
delete [] sendData_Xz; delete [] sendData_XZ; delete [] sendData_yz;
delete [] sendData_yZ; delete [] sendData_Yz; delete [] sendData_YZ;
// Free recvData
delete [] recvData_x; delete [] recvData_y; delete [] recvData_z;
delete [] recvData_X; delete [] recvData_Y; delete [] recvData_Z;
delete [] recvData_xy; delete [] recvData_xY; delete [] recvData_Xy;
delete [] recvData_XY; delete [] recvData_xz; delete [] recvData_xZ;
delete [] recvData_Xz; delete [] recvData_XZ; delete [] recvData_yz;
delete [] recvData_yZ; delete [] recvData_Yz; delete [] recvData_YZ;
// Free id
delete [] id;
// Free the communicator
if ( Comm != MPI_COMM_WORLD && Comm != MPI_COMM_NULL ) {
MPI_Comm_free(&Comm);
}
}
void Domain::initialize( std::shared_ptr<Database> db )
{
d_db = db;
@@ -187,58 +240,337 @@ void Domain::initialize( std::shared_ptr<Database> db )
MPI_Comm_size( Comm, &nprocs );
INSIST(nprocs == nproc[0]*nproc[1]*nproc[2],"Fatal error in processor count!");
}
Domain::~Domain()
{
// Free sendList
delete [] sendList_x; delete [] sendList_y; delete [] sendList_z;
delete [] sendList_X; delete [] sendList_Y; delete [] sendList_Z;
delete [] sendList_xy; delete [] sendList_yz; delete [] sendList_xz;
delete [] sendList_Xy; delete [] sendList_Yz; delete [] sendList_xZ;
delete [] sendList_xY; delete [] sendList_yZ; delete [] sendList_Xz;
delete [] sendList_XY; delete [] sendList_YZ; delete [] sendList_XZ;
// Free sendBuf
delete [] sendBuf_x; delete [] sendBuf_y; delete [] sendBuf_z;
delete [] sendBuf_X; delete [] sendBuf_Y; delete [] sendBuf_Z;
delete [] sendBuf_xy; delete [] sendBuf_yz; delete [] sendBuf_xz;
delete [] sendBuf_Xy; delete [] sendBuf_Yz; delete [] sendBuf_xZ;
delete [] sendBuf_xY; delete [] sendBuf_yZ; delete [] sendBuf_Xz;
delete [] sendBuf_XY; delete [] sendBuf_YZ; delete [] sendBuf_XZ;
// Free recvList
delete [] recvList_x; delete [] recvList_y; delete [] recvList_z;
delete [] recvList_X; delete [] recvList_Y; delete [] recvList_Z;
delete [] recvList_xy; delete [] recvList_yz; delete [] recvList_xz;
delete [] recvList_Xy; delete [] recvList_Yz; delete [] recvList_xZ;
delete [] recvList_xY; delete [] recvList_yZ; delete [] recvList_Xz;
delete [] recvList_XY; delete [] recvList_YZ; delete [] recvList_XZ;
// Free recvBuf
delete [] recvBuf_x; delete [] recvBuf_y; delete [] recvBuf_z;
delete [] recvBuf_X; delete [] recvBuf_Y; delete [] recvBuf_Z;
delete [] recvBuf_xy; delete [] recvBuf_yz; delete [] recvBuf_xz;
delete [] recvBuf_Xy; delete [] recvBuf_Yz; delete [] recvBuf_xZ;
delete [] recvBuf_xY; delete [] recvBuf_yZ; delete [] recvBuf_Xz;
delete [] recvBuf_XY; delete [] recvBuf_YZ; delete [] recvBuf_XZ;
// Free sendData
delete [] sendData_x; delete [] sendData_y; delete [] sendData_z;
delete [] sendData_X; delete [] sendData_Y; delete [] sendData_Z;
delete [] sendData_xy; delete [] sendData_xY; delete [] sendData_Xy;
delete [] sendData_XY; delete [] sendData_xz; delete [] sendData_xZ;
delete [] sendData_Xz; delete [] sendData_XZ; delete [] sendData_yz;
delete [] sendData_yZ; delete [] sendData_Yz; delete [] sendData_YZ;
// Free recvData
delete [] recvData_x; delete [] recvData_y; delete [] recvData_z;
delete [] recvData_X; delete [] recvData_Y; delete [] recvData_Z;
delete [] recvData_xy; delete [] recvData_xY; delete [] recvData_Xy;
delete [] recvData_XY; delete [] recvData_xz; delete [] recvData_xZ;
delete [] recvData_Xz; delete [] recvData_XZ; delete [] recvData_yz;
delete [] recvData_yZ; delete [] recvData_Yz; delete [] recvData_YZ;
// Free id
delete [] id;
// Free the communicator
if ( Comm != MPI_COMM_WORLD && Comm != MPI_COMM_NULL ) {
MPI_Comm_free(&Comm);
}
}
void Domain::Decomp(std::shared_ptr<Database> domain_db )
{
//.......................................................................
// Reading the domain information file
//.......................................................................
int rank_offset = 0;
int RANK = rank();
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz;
int64_t global_Nx,global_Ny,global_Nz;
int64_t i,j,k,n;
int BC=0;
int64_t xStart,yStart,zStart;
int checkerSize;
//int inlet_layers_x, inlet_layers_y, inlet_layers_z;
//int outlet_layers_x, outlet_layers_y, outlet_layers_z;
xStart=yStart=zStart=0;
inlet_layers_x = 0;
inlet_layers_y = 0;
inlet_layers_z = 0;
outlet_layers_x = 0;
outlet_layers_y = 0;
outlet_layers_z = 0;
checkerSize = 32;
// Read domain parameters
auto Filename = domain_db->getScalar<std::string>( "Filename" );
//auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto SIZE = domain_db->getVector<int>( "N" );
auto nproc = domain_db->getVector<int>( "nproc" );
if (domain_db->keyExists( "offset" )){
auto offset = domain_db->getVector<int>( "offset" );
xStart = offset[0];
yStart = offset[1];
zStart = offset[2];
}
if (domain_db->keyExists( "InletLayers" )){
auto InletCount = domain_db->getVector<int>( "InletLayers" );
inlet_layers_x = InletCount[0];
inlet_layers_y = InletCount[1];
inlet_layers_z = InletCount[2];
}
if (domain_db->keyExists( "OutletLayers" )){
auto OutletCount = domain_db->getVector<int>( "OutletLayers" );
outlet_layers_x = OutletCount[0];
outlet_layers_y = OutletCount[1];
outlet_layers_z = OutletCount[2];
}
if (domain_db->keyExists( "checkerSize" )){
checkerSize = domain_db->getScalar<int>( "checkerSize" );
}
else {
checkerSize = SIZE[0];
}
auto ReadValues = domain_db->getVector<int>( "ReadValues" );
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
auto ReadType = domain_db->getScalar<std::string>( "ReadType" );
if (ReadType == "8bit"){
}
else if (ReadType == "16bit"){
}
else{
//printf("INPUT ERROR: Valid ReadType are 8bit, 16bit \n");
ReadType = "8bit";
}
nx = size[0];
ny = size[1];
nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
global_Nx = SIZE[0];
global_Ny = SIZE[1];
global_Nz = SIZE[2];
nprocs=nprocx*nprocy*nprocz;
char *SegData = NULL;
if (RANK==0){
printf("Input media: %s\n",Filename.c_str());
printf("Relabeling %lu values\n",ReadValues.size());
for (int idx=0; idx<ReadValues.size(); idx++){
int oldvalue=ReadValues[idx];
int newvalue=WriteValues[idx];
printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue);
}
// Rank=0 reads the entire segmented data and distributes to worker processes
printf("Dimensions of segmented image: %ld x %ld x %ld \n",global_Nx,global_Ny,global_Nz);
int64_t SIZE = global_Nx*global_Ny*global_Nz;
SegData = new char[SIZE];
if (ReadType == "8bit"){
printf("Reading 8-bit input data \n");
FILE *SEGDAT = fopen(Filename.c_str(),"rb");
if (SEGDAT==NULL) ERROR("Domain.cpp: Error reading segmented data");
size_t ReadSeg;
ReadSeg=fread(SegData,1,SIZE,SEGDAT);
if (ReadSeg != size_t(SIZE)) printf("Domain.cpp: Error reading segmented data \n");
fclose(SEGDAT);
}
else if (ReadType == "16bit"){
printf("Reading 16-bit input data \n");
short int *InputData;
InputData = new short int[SIZE];
FILE *SEGDAT = fopen(Filename.c_str(),"rb");
if (SEGDAT==NULL) ERROR("Domain.cpp: Error reading segmented data");
size_t ReadSeg;
ReadSeg=fread(InputData,2,SIZE,SEGDAT);
if (ReadSeg != size_t(SIZE)) printf("Domain.cpp: Error reading segmented data \n");
fclose(SEGDAT);
for (int n=0; n<SIZE; n++){
SegData[n] = char(InputData[n]);
}
}
printf("Read segmented data from %s \n",Filename.c_str());
if (inlet_layers_x > 0){
// use checkerboard pattern
printf("Checkerboard pattern at x inlet for %i layers \n",inlet_layers_x);
for (int k = 0; k<global_Nz; k++){
for (int j = 0; j<global_Ny; j++){
for (int i = xStart; i < xStart+inlet_layers_x; i++){
if ( (j/checkerSize + k/checkerSize)%2 == 0){
// void checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
}
else{
// solid checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 0;
}
}
}
}
}
if (inlet_layers_y > 0){
printf("Checkerboard pattern at y inlet for %i layers \n",inlet_layers_y);
// use checkerboard pattern
for (int k = 0; k<global_Nz; k++){
for (int j = yStart; i < yStart+inlet_layers_y; j++){
for (int i = 0; i<global_Nx; i++){
if ( (i/checkerSize + k/checkerSize)%2 == 0){
// void checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
}
else{
// solid checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 0;
}
}
}
}
}
if (inlet_layers_z > 0){
printf("Checkerboard pattern at z inlet for %i layers \n",inlet_layers_z);
// use checkerboard pattern
for (int k = zStart; k < zStart+inlet_layers_z; k++){
for (int j = 0; j<global_Ny; j++){
for (int i = 0; i<global_Nx; i++){
if ( (i/checkerSize+j/checkerSize)%2 == 0){
// void checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
}
else{
// solid checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 0;
}
}
}
}
}
if (outlet_layers_x > 0){
// use checkerboard pattern
printf("Checkerboard pattern at x outlet for %i layers \n",outlet_layers_x);
for (int k = 0; k<global_Nz; k++){
for (int j = 0; j<global_Ny; j++){
for (int i = xStart + nx*nprocx - outlet_layers_x; i < xStart + nx*nprocx; i++){
if ( (j/checkerSize + k/checkerSize)%2 == 0){
// void checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
}
else{
// solid checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 0;
}
}
}
}
}
if (outlet_layers_y > 0){
printf("Checkerboard pattern at y outlet for %i layers \n",outlet_layers_y);
// use checkerboard pattern
for (int k = 0; k<global_Nz; k++){
for (int j = yStart + ny*nprocy - outlet_layers_y; i < yStart + ny*nprocy; j++){
for (int i = 0; i<global_Nx; i++){
if ( (i/checkerSize + k/checkerSize)%2 == 0){
// void checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
}
else{
// solid checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 0;
}
}
}
}
}
if (outlet_layers_z > 0){
printf("Checkerboard pattern at z outlet for %i layers \n",outlet_layers_z);
// use checkerboard pattern
for (int k = zStart + nz*nprocz - outlet_layers_z; k < zStart + nz*nprocz; k++){
for (int j = 0; j<global_Ny; j++){
for (int i = 0; i<global_Nx; i++){
if ( (i/checkerSize+j/checkerSize)%2 == 0){
// void checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
}
else{
// solid checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 0;
}
}
}
}
}
}
// Get the rank info
int64_t N = (nx+2)*(ny+2)*(nz+2);
// number of sites to use for periodic boundary condition transition zone
int64_t z_transition_size = (nprocz*nz - (global_Nz - zStart))/2;
if (z_transition_size < 0) z_transition_size=0;
char LocalRankFilename[40];
char *loc_id;
loc_id = new char [(nx+2)*(ny+2)*(nz+2)];
std::vector<int> LabelCount(ReadValues.size(),0);
// Set up the sub-domains
if (RANK==0){
printf("Distributing subdomains across %i processors \n",nprocs);
printf("Process grid: %i x %i x %i \n",nprocx,nprocy,nprocz);
printf("Subdomain size: %i x %i x %i \n",nx,ny,nz);
printf("Size of transition region: %ld \n", z_transition_size);
for (int kp=0; kp<nprocz; kp++){
for (int jp=0; jp<nprocy; jp++){
for (int ip=0; ip<nprocx; ip++){
// rank of the process that gets this subdomain
int rnk = kp*nprocx*nprocy + jp*nprocx + ip;
// Pack and send the subdomain for rnk
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int64_t x = xStart + ip*nx + i-1;
int64_t y = yStart + jp*ny + j-1;
// int64_t z = zStart + kp*nz + k-1;
int64_t z = zStart + kp*nz + k-1 - z_transition_size;
if (x<xStart) x=xStart;
if (!(x<global_Nx)) x=global_Nx-1;
if (y<yStart) y=yStart;
if (!(y<global_Ny)) y=global_Ny-1;
if (z<zStart) z=zStart;
if (!(z<global_Nz)) z=global_Nz-1;
int64_t nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
int64_t nglobal = z*global_Nx*global_Ny+y*global_Nx+x;
loc_id[nlocal] = SegData[nglobal];
}
}
}
// relabel the data
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
n = k*(nx+2)*(ny+2) + j*(nx+2) + i;;
char locval = loc_id[n];
for (int idx=0; idx<ReadValues.size(); idx++){
signed char oldvalue=ReadValues[idx];
signed char newvalue=WriteValues[idx];
if (locval == oldvalue){
loc_id[n] = newvalue;
LabelCount[idx]++;
idx = ReadValues.size();
}
}
}
}
}
if (rnk==0){
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
id[nlocal] = loc_id[nlocal];
}
}
}
}
else{
//printf("Sending data to process %i \n", rnk);
MPI_Send(loc_id,N,MPI_CHAR,rnk,15,Comm);
}
// Write the data for this rank data
sprintf(LocalRankFilename,"ID.%05i",rnk+rank_offset);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(loc_id,1,(nx+2)*(ny+2)*(nz+2),ID);
fclose(ID);
}
}
}
for (int idx=0; idx<ReadValues.size(); idx++){
int label=ReadValues[idx];
int count=LabelCount[idx];
printf("Label=%d, Count=%d \n",label,count);
}
}
else{
// Recieve the subdomain from rank = 0
//printf("Ready to recieve data %i at process %i \n", N,rank);
MPI_Recv(id,N,MPI_CHAR,0,15,Comm,MPI_STATUS_IGNORE);
}
MPI_Barrier(Comm);
}
/********************************************************
* Initialize communication *

1
common/Domain.h
View File

@@ -190,6 +190,7 @@ public: // Public variables (need to create accessors instead)
signed char *id;
void ReadIDs();
void Decomp(std::shared_ptr<Database> domain_db);
void CommunicateMeshHalo(DoubleArray &Mesh);
void CommInit();
int PoreCount();

96
models/ColorModel.cpp
View File

@@ -75,7 +75,7 @@ void ScaLBL_ColorModel::ReadParams(string filename){
// read the input database
db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" );
color_db = db->getDatabase( "Color" );
color_db = db->getDatabase( "Color" );
analysis_db = db->getDatabase( "Analysis" );
// set defaults
@@ -166,8 +166,14 @@ void ScaLBL_ColorModel::SetDomain(){
}
void ScaLBL_ColorModel::ReadInput(){
size_t readID;
Mask->ReadIDs();
if (domain_db->keyExists( "Filename" )){
Mask->Decomp(domain_db);
}
else{
size_t readID;
Mask->ReadIDs();
}
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
sprintf(LocalRankString,"%05d",rank);
@@ -648,6 +654,7 @@ void ScaLBL_ColorModel::Run(){
double volA = Averages->gnb.V;
volA /= Dm->Volume;
volB /= Dm->Volume;;
initial_volume = volA*Dm->Volume;
double vA_x = Averages->gnb.Px/Averages->gnb.M;
double vA_y = Averages->gnb.Py/Averages->gnb.M;
double vA_z = Averages->gnb.Pz/Averages->gnb.M;
@@ -681,10 +688,9 @@ void ScaLBL_ColorModel::Run(){
isSteady = true;
if ( isSteady ){
initial_volume = volA*Dm->Volume;
MORPH_ADAPT = true;
CURRENT_MORPH_TIMESTEPS=0;
delta_volume_target = Dm->Volume*morph_delta; //*volA ???? // set target volume change
delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
Averages->Full();
Averages->Write(timestep);
analysis.WriteVisData( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
@@ -763,8 +769,7 @@ void ScaLBL_ColorModel::Run(){
else if (USE_SEED){
delta_volume = volA*Dm->Volume - initial_volume;
CURRENT_MORPH_TIMESTEPS += analysis_interval;
double effective_seed_water = seed_water*(1.0+volB/volA);
double massChange = SeedPhaseField(effective_seed_water);
double massChange = SeedPhaseField(seed_water);
if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target);
}
else if (USE_MORPHOPEN_OIL){
@@ -984,7 +989,6 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
double mass_loss =0.f;
double count =0.f;
double *Aq_tmp, *Bq_tmp;
double SEED_THRESHOLD = -0.9;
Aq_tmp = new double [7*Np];
Bq_tmp = new double [7*Np];
@@ -1012,58 +1016,54 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
}
}
*/
double oil_value = 0.0;
double water_value = 1.0;
for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
double phase_id = (dA - dB) / (dA + dB);
double random_value = double(rand())/ RAND_MAX;
if (phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){
Aq_tmp[n] = 0.3333333333333333*oil_value;
Aq_tmp[n+Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value;
if (phase_id > 0.0){
Aq_tmp[n] -= 0.3333333333333333*random_value;
Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
Bq_tmp[n] = 0.3333333333333333*water_value;
Bq_tmp[n+Np] = 0.1111111111111111*water_value;
Bq_tmp[n+2*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+3*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+4*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+5*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+6*Np] = 0.1111111111111111*water_value;
mass_loss += oil_value - dA;
count++;
Bq_tmp[n] += 0.3333333333333333*random_value;
Bq_tmp[n+Np] += 0.1111111111111111*random_value;
Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
}
mass_loss += random_value*seed_water_in_oil;
}
for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
double phase_id = (dA - dB) / (dA + dB);
double random_value = double(rand())/ RAND_MAX;
if (phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){
Aq_tmp[n] = 0.3333333333333333*oil_value;
Aq_tmp[n+Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value;
Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value;
if (phase_id > 0.0){
Aq_tmp[n] -= 0.3333333333333333*random_value;
Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
Bq_tmp[n] = 0.3333333333333333*water_value;
Bq_tmp[n+Np] = 0.1111111111111111*water_value;
Bq_tmp[n+2*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+3*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+4*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+5*Np] = 0.1111111111111111*water_value;
Bq_tmp[n+6*Np] = 0.1111111111111111*water_value;
mass_loss += oil_value - dA;
count++;
Bq_tmp[n] += 0.3333333333333333*random_value;
Bq_tmp[n+Np] += 0.1111111111111111*random_value;
Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
}
mass_loss += random_value*seed_water_in_oil;
}
count= sumReduce( Dm->Comm, count);

30
sample_scripts/configure_huckleberry Executable file
View File

@@ -0,0 +1,30 @@
# load the module for cmake
module load cmake
module load gcc/7.3.0 openmpi/3.1.2 cuda
module load szip
module load zlib
module list
module load hdf5/1.8.12
module load silo
cmake \
-D CMAKE_C_COMPILER:PATH=mpicc \
-D CMAKE_CXX_COMPILER:PATH=mpicxx \
-D CMAKE_C_FLAGS="-fPIC" \
-D CMAKE_CXX_FLAGS="-fPIC" \
-D MPI_COMPILER:BOOL=TRUE \
-D MPIEXEC=mpirun \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D CMAKE_BUILD_TYPE:STRING=Release \
-D USE_CUDA=1 \
-D CMAKE_CUDA_FLAGS="-arch sm_60" \
-D CMAKE_CUDA_HOST_COMPILER="$GCC_BIN/gcc" \
-D USE_HDF5=1 \
-D HDF5_DIRECTORY=$HDF5_DIR \
-D USE_SILO=1 \
-D SILO_DIRECTORY=$SILO_DIR \
-D USE_TIMER=0 \
~/LBPM-WIA
make VERBOSE=1 -j1 && make install

2
tests/CMakeLists.txt
View File

@@ -39,6 +39,7 @@ CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_
# Add the tests
ADD_LBPM_TEST( pmmc_cylinder )
ADD_LBPM_TEST( TestSubphase )
ADD_LBPM_TEST( TestTorus )
ADD_LBPM_TEST( TestTorusEvolve )
ADD_LBPM_TEST( TestTopo3D )
@@ -67,6 +68,7 @@ ADD_LBPM_TEST_PARALLEL( TestSegDist 8 )
ADD_LBPM_TEST_PARALLEL( TestCommD3Q19 8 )
ADD_LBPM_TEST_1_2_4( testCommunication )
ADD_LBPM_TEST( TestWriter )
ADD_LBPM_TEST( TestDatabase )
IF ( USE_NETCDF )
ADD_LBPM_TEST_PARALLEL( TestNetcdf 8 )
ADD_LBPM_EXECUTABLE( lbpm_uCT_pp )

71
tests/TestDatabase.cpp Normal file
View File

@@ -0,0 +1,71 @@
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <exception>
#include <stdexcept>
#include <fstream>
#include <iostream>
#include <memory>
#include "common/UnitTest.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "common/Database.h"
#include "ProfilerApp.h"
// Main
int main(int argc, char **argv)
{
int rank,nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
Utilities::setAbortBehavior(true,2);
Utilities::setErrorHandlers();
UnitTest ut;
int err=0;
int BC=2;
int npx=1; int npy=2; int npz=4;
int nx=32; int ny=34; int nz=35;
std::vector<std::string> List;
List.push_back("name1");
List.push_back("name2");
// write a simple database and test that it can be read by LBPM
auto db = std::make_shared<Database>( );
db->putScalar<int>( "BC", BC );
db->putVector<int>( "nproc", { npx, npy, npz } );
db->putVector<int>( "n", { nx, ny, nz } );
db->putVector<std::string>( "Files", List);
std::ofstream OutStream("test.db");
// db->putDatabase();
OutStream << "Domain { \n";
db->print(OutStream, " ");
OutStream << "} \n";
printf("TestDatbase: writing test file\n");
OutStream.close();
auto new_db = std::make_shared<Database>( "test.db" );
auto domain_db = new_db->getDatabase( "Domain" );
if (domain_db->keyExists( "BC" )){
auto newBC = domain_db->getScalar<int>( "BC" );
if (newBC != BC){
err=1;
printf("TestDatbase: getScalar failed! \n");
}
}
if (err==0) printf("TestDatabase: succeeded!\n");
else printf("TestDatabase: errors encountered \n");
// Finished
PROFILE_SAVE("TestDatabase",true);
MPI_Barrier(comm);
MPI_Finalize();
return err;
}

146
tests/TestSubphase.cpp Normal file
View File

@@ -0,0 +1,146 @@
// Sequential blob analysis
// Reads parallel simulation data and performs connectivity analysis
// and averaging on a blob-by-blob basis
// James E. McClure 2014
#include <iostream>
#include <math.h>
#include "common/Communication.h"
#include "analysis/analysis.h"
#include "analysis/SubPhase.h"
std::shared_ptr<Database> loadInputs( int nprocs )
{
//auto db = std::make_shared<Database>( "Domain.in" );
auto db = std::make_shared<Database>();
db->putScalar<int>( "BC", 0 );
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 100, 100, 100 } );
db->putScalar<int>( "nspheres", 1 );
db->putVector<double>( "L", { 1, 1, 1 } );
return db;
}
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);
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
if ( rank==0 ) {
printf("-----------------------------------------------------------\n");
printf("Unit test for torus (Euler-Poincarie characteristic) \n");
printf("-----------------------------------------------------------\n");
}
//.......................................................................
// Reading the domain information file
//.......................................................................
int i,j,k,n;
// Load inputs
auto db = loadInputs( nprocs );
int Nx = db->getVector<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "nproc" )[2];
if (rank==0){
printf("********************************************************\n");
printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
printf("********************************************************\n");
}
// Get the rank info
std::shared_ptr<Domain> Dm(new Domain(db,comm));
// const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
std::shared_ptr<SubPhase> Averages(new SubPhase(Dm));
Nx += 2;
Ny += 2;
Nz += 2;
int N = Nx*Ny*Nz;
//.......................................................................
for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
Dm->id[n] = 1;
}
}
}
//.......................................................................
Dm->CommInit(); // Initialize communications for domains
//.......................................................................
//.......................................................................
// Assign the phase ID field based and the signed distance
//.......................................................................
double R1,R2;
double CX,CY,CZ; //CY1,CY2;
CX=Nx*nprocx*0.5;
CY=Ny*nprocy*0.5;
CZ=Nz*nprocz*0.5;
R1 = Nx*nprocx*0.2; // middle radius
R2 = Nx*nprocx*0.1; // donut thickness
//
//CY1=Nx*nprocx*0.5+R1;
//CY2=Ny*nprocy*0.5-R1;
double x,y,z;
if (rank==0) printf("Initializing the system \n");
for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
// global position relative to center
x = Dm->iproc()*Nx+i - CX;
y = Dm->jproc()*Ny+j - CY;
z = Dm->kproc()*Nz+k - CZ;
// Shrink the sphere sizes by two voxels to make sure they don't touch
Averages->SDs(i,j,k) = 100.0;
//..............................................................................
// Single torus
Averages->Phi(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2;
// Double torus
/* y = Dm->jproc()*Ny+j - CY1;
//z = Dm->kproc()*Nz+k - CZ +R1;
Averages->Phi(i,j,k) = sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2;
y = Dm->jproc()*Ny+j - CY2;
//z = Dm->kproc()*Nz+k - CZ-R1;
Averages->Phi(i,j,k) = min(Averages->Phi(i,j,k),
sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z) - R2);
*///..............................................................................
//Averages->Phi(i,j,k) = - Averages->Phi(i,j,k);
if (Averages->Phi(i,j,k) > 0.0){
Dm->id[n] = 2;
}
else{
Dm->id[n] = 1;
}
}
}
}
if (rank==0) printf("Aggregating labels \n");
Averages->AggregateLabels("phase_id.raw");
// Averages->Reduce();
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Barrier(comm);
MPI_Finalize();
return 0;
}

67
tests/lbpm_serial_decomp.cpp
View File

@@ -55,12 +55,16 @@ int main(int argc, char **argv)
int64_t xStart,yStart,zStart;
int checkerSize;
int inlet_count_x, inlet_count_y, inlet_count_z;
int outlet_count_x, outlet_count_y, outlet_count_z;
// char fluidValue,solidValue;
xStart=yStart=zStart=0;
inlet_count_x = 0;
inlet_count_y = 0;
inlet_count_z = 0;
outlet_count_x = 0;
outlet_count_y = 0;
outlet_count_z = 0;
checkerSize = 32;
// read the input database
auto db = std::make_shared<Database>( filename );
@@ -84,6 +88,12 @@ int main(int argc, char **argv)
inlet_count_y = InletCount[1];
inlet_count_z = InletCount[2];
}
if (domain_db->keyExists( "OutletLayers" )){
auto OutletCount = domain_db->getVector<int>( "OutletLayers" );
outlet_count_x = OutletCount[0];
outlet_count_y = OutletCount[1];
outlet_count_z = OutletCount[2];
}
if (domain_db->keyExists( "checkerSize" )){
checkerSize = domain_db->getScalar<int>( "checkerSize" );
}
@@ -211,6 +221,63 @@ int main(int argc, char **argv)
}
}
if (outlet_count_x > 0){
// use checkerboard pattern
printf("Checkerboard pattern at x outlet for %i layers \n",outlet_count_x);
for (int k = 0; k<Nz; k++){
for (int j = 0; j<Ny; j++){
for (int i = xStart + nx*nprocx - outlet_count_x; i < xStart + nx*nprocx; i++){
if ( (j/checkerSize + k/checkerSize)%2 == 0){
// void checkers
SegData[k*Nx*Ny+j*Nx+i] = 2;
}
else{
// solid checkers
SegData[k*Nx*Ny+j*Nx+i] = 0;
}
}
}
}
}
if (outlet_count_y > 0){
printf("Checkerboard pattern at y outlet for %i layers \n",outlet_count_y);
// use checkerboard pattern
for (int k = 0; k<Nz; k++){
for (int j = yStart + ny*nprocy - outlet_count_y; i < yStart + ny*nprocy; j++){
for (int i = 0; i<Nx; i++){
if ( (i/checkerSize + k/checkerSize)%2 == 0){
// void checkers
SegData[k*Nx*Ny+j*Nx+i] = 2;
}
else{
// solid checkers
SegData[k*Nx*Ny+j*Nx+i] = 0;
}
}
}
}
}
if (outlet_count_z > 0){
printf("Checkerboard pattern at z outlet for %i layers \n",outlet_count_z);
// use checkerboard pattern
for (int k = zStart + nz*nprocz - outlet_count_z; k < zStart + nz*nprocz; k++){
for (int j = 0; j<Ny; j++){
for (int i = 0; i<Nx; i++){
if ( (i/checkerSize+j/checkerSize)%2 == 0){
// void checkers
SegData[k*Nx*Ny+j*Nx+i] = 2;
}
else{
// solid checkers
SegData[k*Nx*Ny+j*Nx+i] = 0;
}
}
}
}
}
// Get the rank info
int64_t N = (nx+2)*(ny+2)*(nz+2);