OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
57a723b11e886e8a0e98e97b04d5b75d8a8de359
LBPM/common/Domain.cpp
Rex Zhe Li 57a723b11e fix a minor data type bug
2021-06-18 00:13:17 -04:00

1457 lines
54 KiB
C++
Raw Blame History

// Created by James McClure
// Copyright 2008-2020
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include <time.h>
#include <exception> // std::exception
#include <stdexcept>
#include "common/Domain.h"
#include "common/Array.h"
#include "common/Utilities.h"
#include "common/MPI.h"
#include "common/Communication.h"
// Inline function to read line without a return argument
static inline void fgetl( char * str, int num, FILE * stream )
{
char* ptr = fgets( str, num, stream );
if ( 0 ) {char *temp = (char *)&ptr; temp++;}
}
/********************************************************
* Constructors *
********************************************************/
Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
double lx, double ly, double lz, int BC):
database(nullptr), Nx(0), Ny(0), Nz(0),
Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0), voxel_length(1),
Comm( Utilities::MPI( MPI_COMM_WORLD).dup() ),
inlet_layers_x(0), inlet_layers_y(0), inlet_layers_z(0),
inlet_layers_phase(1),outlet_layers_phase(2)
{
NULL_USE( rnk );
NULL_USE( npy );
NULL_USE( npz );
// set up the neighbor ranks
int myrank = Comm.getRank();
rank_info = RankInfoStruct( myrank, rank_info.nx, rank_info.ny, rank_info.nz );
Comm.barrier();
auto db = std::make_shared<Database>( );
db->putScalar<int>( "BC", BC );
db->putVector<int>( "nproc", { npx, npx, npx } );
db->putVector<int>( "n", { nx, ny, nz } );
db->putScalar<int>( "nspheres", 0 );
db->putVector<double>( "L", { lx, ly, lz } );
initialize( db );
}
Domain::Domain( std::shared_ptr<Database> db, const Utilities::MPI& Communicator):
database(db), Nx(0), Ny(0), Nz(0),
Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(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),
inlet_layers_phase(1),outlet_layers_phase(2)
{
Comm = Communicator.dup();
// set up the neighbor ranks
int myrank = Comm.getRank();
initialize( db );
rank_info = RankInfoStruct( myrank, rank_info.nx, rank_info.ny, rank_info.nz );
Comm.barrier();
}
/********************************************************
* Destructor *
********************************************************/
Domain::~Domain()
{
}
/********************************************************
* Initialization *
********************************************************/
void Domain::initialize( std::shared_ptr<Database> db )
{
d_db = db;
auto nproc = d_db->getVector<int>("nproc");
auto n = d_db->getVector<int>("n");
ASSERT( n.size() == 3u );
ASSERT( nproc.size() == 3u );
int nx = n[0];
int ny = n[1];
int nz = n[2];
if (d_db->keyExists( "InletLayers" )){
auto InletCount = d_db->getVector<int>( "InletLayers" );
inlet_layers_x = InletCount[0];
inlet_layers_y = InletCount[1];
inlet_layers_z = InletCount[2];
}
if (d_db->keyExists( "OutletLayers" )){
auto OutletCount = d_db->getVector<int>( "OutletLayers" );
outlet_layers_x = OutletCount[0];
outlet_layers_y = OutletCount[1];
outlet_layers_z = OutletCount[2];
}
if (d_db->keyExists( "InletLayersPhase" )){
inlet_layers_phase = d_db->getScalar<int>( "InletLayersPhase" );
}
if (d_db->keyExists( "OutletLayersPhase" )){
outlet_layers_phase = d_db->getScalar<int>( "OutletLayersPhase" );
}
voxel_length = 1.0;
if (d_db->keyExists( "voxel_length" )){
voxel_length = d_db->getScalar<double>("voxel_length");
}
else if (d_db->keyExists( "L" )){
auto Length = d_db->getVector<double>("L");
Lx = Length[0];
Ly = Length[1];
Lz = Length[2];
voxel_length = Lx/(nx*nproc[0]);
}
Lx = nx*nproc[0]*voxel_length;
Ly = ny*nproc[1]*voxel_length;
Lz = nz*nproc[2]*voxel_length;
Nx = nx+2;
Ny = ny+2;
Nz = nz+2;
// Initialize ranks
int myrank = Comm.getRank();
rank_info = RankInfoStruct(myrank,nproc[0],nproc[1],nproc[2]);
// inlet layers only apply to lower part of domain
if (rank_info.ix > 0) inlet_layers_x = 0;
if (rank_info.jy > 0) inlet_layers_y = 0;
if (rank_info.kz > 0) inlet_layers_z = 0;
// outlet layers only apply to top part of domain
if (rank_info.ix < nproc[0]-1 ) outlet_layers_x = 0;
if (rank_info.jy < nproc[1]-1) outlet_layers_y = 0;
if (rank_info.kz < nproc[2]-1) outlet_layers_z = 0;
// Fill remaining variables
N = Nx*Ny*Nz;
Volume = nx*ny*nz*nproc[0]*nproc[1]*nproc[2]*1.0;
if (myrank==0) printf("voxel length = %f micron \n", voxel_length);
id = std::vector<signed char>( N, 0 );
BoundaryCondition = d_db->getScalar<int>("BC");
int nprocs = Comm.getSize();
INSIST(nprocs == nproc[0]*nproc[1]*nproc[2],"Fatal error in processor count!");
}
/********************************************************
* Get send/recv lists *
********************************************************/
const std::vector<int>& Domain::getRecvList( const char* dir ) const
{
if ( dir[0] == 'x' ) {
if ( dir[1] == 0 )
return recvList_x;
else if ( dir[1] == 'y' )
return recvList_xy;
else if ( dir[1] == 'Y' )
return recvList_xY;
else if ( dir[1] == 'z' )
return recvList_xz;
else if ( dir[1] == 'Z' )
return recvList_xZ;
} else if ( dir[0] == 'y' ) {
if ( dir[1] == 0 )
return recvList_y;
else if ( dir[1] == 'z' )
return recvList_yz;
else if ( dir[1] == 'Z' )
return recvList_yZ;
} else if ( dir[0] == 'z' ) {
if ( dir[1] == 0 )
return recvList_z;
} else if ( dir[0] == 'X' ) {
if ( dir[1] == 0 )
return recvList_X;
else if ( dir[1] == 'y' )
return recvList_Xy;
else if ( dir[1] == 'Y' )
return recvList_XY;
else if ( dir[1] == 'z' )
return recvList_Xz;
else if ( dir[1] == 'Z' )
return recvList_XZ;
} else if ( dir[0] == 'Y' ) {
if ( dir[1] == 0 )
return recvList_Y;
else if ( dir[1] == 'z' )
return recvList_Yz;
else if ( dir[1] == 'Z' )
return recvList_YZ;
} else if ( dir[0] == 'Z' ) {
if ( dir[1] == 0 )
return recvList_Z;
}
throw std::logic_error("Internal error");
}
const std::vector<int>& Domain::getSendList( const char* dir ) const
{
if ( dir[0] == 'x' ) {
if ( dir[1] == 0 )
return sendList_x;
else if ( dir[1] == 'y' )
return sendList_xy;
else if ( dir[1] == 'Y' )
return sendList_xY;
else if ( dir[1] == 'z' )
return sendList_xz;
else if ( dir[1] == 'Z' )
return sendList_xZ;
} else if ( dir[0] == 'y' ) {
if ( dir[1] == 0 )
return sendList_y;
else if ( dir[1] == 'z' )
return sendList_yz;
else if ( dir[1] == 'Z' )
return sendList_yZ;
} else if ( dir[0] == 'z' ) {
if ( dir[1] == 0 )
return sendList_z;
} else if ( dir[0] == 'X' ) {
if ( dir[1] == 0 )
return sendList_X;
else if ( dir[1] == 'y' )
return sendList_Xy;
else if ( dir[1] == 'Y' )
return sendList_XY;
else if ( dir[1] == 'z' )
return sendList_Xz;
else if ( dir[1] == 'Z' )
return sendList_XZ;
} else if ( dir[0] == 'Y' ) {
if ( dir[1] == 0 )
return sendList_Y;
else if ( dir[1] == 'z' )
return sendList_Yz;
else if ( dir[1] == 'Z' )
return sendList_YZ;
} else if ( dir[0] == 'Z' ) {
if ( dir[1] == 0 )
return sendList_Z;
}
throw std::logic_error("Internal error");
}
/********************************************************
* Decomp *
********************************************************/
void Domain::Decomp( const std::string& Filename )
{
//.......................................................................
// 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;
int64_t xStart,yStart,zStart;
int checkerSize;
bool USE_CHECKER = false;
//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;
inlet_layers_phase=1;
outlet_layers_phase=2;
checkerSize = 32;
// Read domain parameters
//auto Filename = database->getScalar<std::string>( "Filename" );
//auto L = database->getVector<double>( "L" );
auto size = database->getVector<int>( "n" );
auto SIZE = database->getVector<int>( "N" );
auto nproc = database->getVector<int>( "nproc" );
if (database->keyExists( "offset" )){
auto offset = database->getVector<int>( "offset" );
xStart = offset[0];
yStart = offset[1];
zStart = offset[2];
}
if (database->keyExists( "InletLayers" )){
auto InletCount = database->getVector<int>( "InletLayers" );
inlet_layers_x = InletCount[0];
inlet_layers_y = InletCount[1];
inlet_layers_z = InletCount[2];
}
if (database->keyExists( "OutletLayers" )){
auto OutletCount = database->getVector<int>( "OutletLayers" );
outlet_layers_x = OutletCount[0];
outlet_layers_y = OutletCount[1];
outlet_layers_z = OutletCount[2];
}
if (database->keyExists( "checkerSize" )){
checkerSize = database->getScalar<int>( "checkerSize" );
USE_CHECKER = true;
}
else {
checkerSize = SIZE[0];
}
if (database->keyExists( "InletLayersPhase" )){
inlet_layers_phase = database->getScalar<int>( "InletLayersPhase" );
}
if (database->keyExists( "OutletLayersPhase" )){
outlet_layers_phase = database->getScalar<int>( "OutletLayersPhase" );
}
auto ReadValues = database->getVector<int>( "ReadValues" );
auto WriteValues = database->getVector<int>( "WriteValues" );
auto ReadType = database->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 (size_t 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());
// relabel the data
std::vector<long int> LabelCount(ReadValues.size(),0);
for (int k = 0; k<global_Nz; k++){
for (int j = 0; j<global_Ny; j++){
for (int i = 0; i<global_Nx; i++){
n = k*global_Nx*global_Ny+j*global_Nx+i;
//char locval = loc_id[n];
signed char locval = SegData[n];
for (size_t idx=0; idx<ReadValues.size(); idx++){
signed char oldvalue=ReadValues[idx];
signed char newvalue=WriteValues[idx];
if (locval == oldvalue){
SegData[n] = newvalue;
LabelCount[idx]++;
idx = ReadValues.size();
}
}
}
}
}
for (size_t idx=0; idx<ReadValues.size(); idx++){
long int label=ReadValues[idx];
long int count=LabelCount[idx];
printf("Label=%ld, Count=%ld \n",label,count);
}
if (USE_CHECKER) {
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; j < 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, saturated with phase label=%i \n",inlet_layers_z,inlet_layers_phase);
// 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;
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = inlet_layers_phase;
}
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; j < 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, saturated with phase label=%i \n",outlet_layers_z,outlet_layers_phase);
// 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;
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = outlet_layers_phase;
}
else{
// solid checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 0;
}
}
}
}
}
}
else {
if (inlet_layers_z > 0){
printf("Mixed reflection pattern at z inlet for %i layers, saturated with phase label=%i \n",inlet_layers_z,inlet_layers_phase);
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++){
signed char local_id = SegData[k*global_Nx*global_Ny+j*global_Nx+i];
signed char reflection_id = SegData[(zStart + nz*nprocz - 1)*global_Nx*global_Ny+j*global_Nx+i];
if ( local_id < 1 && reflection_id > 0){
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = reflection_id;
}
}
}
}
}
if (outlet_layers_z > 0){
printf("Mixed reflection pattern at z outlet for %i layers, saturated with phase label=%i \n",outlet_layers_z,outlet_layers_phase);
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++){
signed char local_id = SegData[k*global_Nx*global_Ny+j*global_Nx+i];
signed char reflection_id = SegData[zStart*global_Nx*global_Ny+j*global_Nx+i];
if ( local_id < 1 && reflection_id > 0){
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = reflection_id;
}
}
}
}
}
}
}
// 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;
// Set up the sub-domains
if (RANK==0){
printf("Distributing subdomains across %i processors \n",nprocs);
printf("Process grid: %i x %i x %i \n",nprocx,nprocy,nprocz);
printf("Subdomain size: %i x %i x %i \n",nx,ny,nz);
printf("Size of transition region: %ld \n", z_transition_size);
auto loc_id = new char [(nx+2)*(ny+2)*(nz+2)];
for (int kp=0; kp<nprocz; kp++){
for (int jp=0; jp<nprocy; jp++){
for (int ip=0; ip<nprocx; ip++){
// 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];
}
}
}
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);
Comm.send(loc_id,N,rnk,15);
}
// Write the data for this rank data
char LocalRankFilename[40];
sprintf(LocalRankFilename,"ID.%05i",rnk+rank_offset);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(loc_id,1,(nx+2)*(ny+2)*(nz+2),ID);
fclose(ID);
}
}
}
delete [] loc_id;
} else {
// Recieve the subdomain from rank = 0
//printf("Ready to recieve data %i at process %i \n", N,rank);
Comm.recv(id.data(),N,0,15);
}
Comm.barrier();
ComputePorosity();
delete [] SegData;
}
void Domain::ComputePorosity(){
// Compute the porosity
double sum;
double sum_local=0.0;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocx()*nprocy()*nprocz());
if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
//.........................................................
for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-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;
if (id[n] > 0){
sum_local+=1.0;
}
}
}
}
sum = Comm.sumReduce(sum_local);
porosity = sum*iVol_global;
if (rank()==0) printf("Media porosity = %f \n",porosity);
//.........................................................
}
void Domain::AggregateLabels( const std::string& filename ){
int nx = Nx;
int ny = Ny;
int nz = Nz;
int npx = nprocx();
int npy = nprocy();
int npz = nprocz();
int ipx = iproc();
int ipy = jproc();
int ipz = kproc();
int nprocs = nprocx()*nprocy()*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);
auto LocalID = new signed char [local_size];
//printf("aggregate labels: local size=%i, global size = %i",local_size, full_size);
// assign the ID for the local sub-region
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 = id[n];
LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val;
}
}
}
Comm.barrier();
// populate the FullID
if (rank() == 0){
auto FullID = new signed char [full_size];
// first handle local ID for rank 0
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
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;
Comm.recv(LocalID,local_size,rnk,tag);
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 = fopen(filename.c_str(),"wb");
fwrite(FullID,1,full_size,OUTFILE);
fclose(OUTFILE);
delete [] FullID;
}
else{
// send LocalID to rank=0
int tag = 15+ rank();
int dstrank = 0;
Comm.send(LocalID,local_size,dstrank,tag);
}
delete [] LocalID;
Comm.barrier();
}
/********************************************************
* Initialize communication *
********************************************************/
void Domain::CommInit()
{
int i,j,k,n;
int sendtag = 21;
int recvtag = 21;
//......................................................................................
int sendCount_x, sendCount_y, sendCount_z, sendCount_X, sendCount_Y, sendCount_Z;
int sendCount_xy, sendCount_yz, sendCount_xz, sendCount_Xy, sendCount_Yz, sendCount_xZ;
int sendCount_xY, sendCount_yZ, sendCount_Xz, sendCount_XY, sendCount_YZ, sendCount_XZ;
sendCount_x = sendCount_y = sendCount_z = sendCount_X = sendCount_Y = sendCount_Z = 0;
sendCount_xy = sendCount_yz = sendCount_xz = sendCount_Xy = sendCount_Yz = sendCount_xZ = 0;
sendCount_xY = sendCount_yZ = sendCount_Xz = sendCount_XY = sendCount_YZ = sendCount_XZ = 0;
//......................................................................................
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
// Check the phase ID
if (id[k*Nx*Ny+j*Nx+i] > 0){
// Counts for the six faces
if (i==1) sendCount_x++;
if (j==1) sendCount_y++;
if (k==1) sendCount_z++;
if (i==Nx-2) sendCount_X++;
if (j==Ny-2) sendCount_Y++;
if (k==Nz-2) sendCount_Z++;
// Counts for the twelve edges
if (i==1 && j==1) sendCount_xy++;
if (i==1 && j==Ny-2) sendCount_xY++;
if (i==Nx-2 && j==1) sendCount_Xy++;
if (i==Nx-2 && j==Ny-2) sendCount_XY++;
if (i==1 && k==1) sendCount_xz++;
if (i==1 && k==Nz-2) sendCount_xZ++;
if (i==Nx-2 && k==1) sendCount_Xz++;
if (i==Nx-2 && k==Nz-2) sendCount_XZ++;
if (j==1 && k==1) sendCount_yz++;
if (j==1 && k==Nz-2) sendCount_yZ++;
if (j==Ny-2 && k==1) sendCount_Yz++;
if (j==Ny-2 && k==Nz-2) sendCount_YZ++;
}
}
}
}
// allocate send lists
sendList_x.resize( sendCount_x, 0 );
sendList_y.resize( sendCount_y, 0 );
sendList_z.resize( sendCount_z, 0 );
sendList_X.resize( sendCount_X, 0 );
sendList_Y.resize( sendCount_Y, 0 );
sendList_Z.resize( sendCount_Z, 0 );
sendList_xy.resize( sendCount_xy, 0 );
sendList_yz.resize( sendCount_yz, 0 );
sendList_xz.resize( sendCount_xz, 0 );
sendList_Xy.resize( sendCount_Xy, 0 );
sendList_Yz.resize( sendCount_Yz, 0 );
sendList_xZ.resize( sendCount_xZ, 0 );
sendList_xY.resize( sendCount_xY, 0 );
sendList_yZ.resize( sendCount_yZ, 0 );
sendList_Xz.resize( sendCount_Xz, 0 );
sendList_XY.resize( sendCount_XY, 0 );
sendList_YZ.resize( sendCount_YZ, 0 );
sendList_XZ.resize( sendCount_XZ, 0 );
// Populate the send list
sendCount_x = sendCount_y = sendCount_z = sendCount_X = sendCount_Y = sendCount_Z = 0;
sendCount_xy = sendCount_yz = sendCount_xz = sendCount_Xy = sendCount_Yz = sendCount_xZ = 0;
sendCount_xY = sendCount_yZ = sendCount_Xz = sendCount_XY = sendCount_YZ = sendCount_XZ = 0;
for (k=1; k<Nz-1; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
// Local value to send
n = k*Nx*Ny+j*Nx+i;
if (id[n] > 0){
// Counts for the six faces
if (i==1) sendList_x[sendCount_x++]=n;
if (j==1) sendList_y[sendCount_y++]=n;
if (k==1) sendList_z[sendCount_z++]=n;
if (i==Nx-2) sendList_X[sendCount_X++]=n;
if (j==Ny-2) sendList_Y[sendCount_Y++]=n;
if (k==Nz-2) sendList_Z[sendCount_Z++]=n;
// Counts for the twelve edges
if (i==1 && j==1) sendList_xy[sendCount_xy++]=n;
if (i==1 && j==Ny-2) sendList_xY[sendCount_xY++]=n;
if (i==Nx-2 && j==1) sendList_Xy[sendCount_Xy++]=n;
if (i==Nx-2 && j==Ny-2) sendList_XY[sendCount_XY++]=n;
if (i==1 && k==1) sendList_xz[sendCount_xz++]=n;
if (i==1 && k==Nz-2) sendList_xZ[sendCount_xZ++]=n;
if (i==Nx-2 && k==1) sendList_Xz[sendCount_Xz++]=n;
if (i==Nx-2 && k==Nz-2) sendList_XZ[sendCount_XZ++]=n;
if (j==1 && k==1) sendList_yz[sendCount_yz++]=n;
if (j==1 && k==Nz-2) sendList_yZ[sendCount_yZ++]=n;
if (j==Ny-2 && k==1) sendList_Yz[sendCount_Yz++]=n;
if (j==Ny-2 && k==Nz-2) sendList_YZ[sendCount_YZ++]=n;
}
}
}
}
//......................................................................................
int recvCount_x, recvCount_y, recvCount_z, recvCount_X, recvCount_Y, recvCount_Z;
int recvCount_xy, recvCount_yz, recvCount_xz, recvCount_Xy, recvCount_Yz, recvCount_xZ;
int recvCount_xY, recvCount_yZ, recvCount_Xz, recvCount_XY, recvCount_YZ, recvCount_XZ;
req1[0] = Comm.Isend(&sendCount_x,1,rank_x(),sendtag+0);
req2[0] = Comm.Irecv(&recvCount_X,1,rank_X(),recvtag+0);
req1[1] = Comm.Isend(&sendCount_X,1,rank_X(),sendtag+1);
req2[1] = Comm.Irecv(&recvCount_x,1,rank_x(),recvtag+1);
req1[2] = Comm.Isend(&sendCount_y,1,rank_y(),sendtag+2);
req2[2] = Comm.Irecv(&recvCount_Y,1,rank_Y(),recvtag+2);
req1[3] = Comm.Isend(&sendCount_Y,1,rank_Y(),sendtag+3);
req2[3] = Comm.Irecv(&recvCount_y,1,rank_y(),recvtag+3);
req1[4] = Comm.Isend(&sendCount_z,1,rank_z(),sendtag+4);
req2[4] = Comm.Irecv(&recvCount_Z,1,rank_Z(),recvtag+4);
req1[5] = Comm.Isend(&sendCount_Z,1,rank_Z(),sendtag+5);
req2[5] = Comm.Irecv(&recvCount_z,1,rank_z(),recvtag+5);
req1[6] = Comm.Isend(&sendCount_xy,1,rank_xy(),sendtag+6);
req2[6] = Comm.Irecv(&recvCount_XY,1,rank_XY(),recvtag+6);
req1[7] = Comm.Isend(&sendCount_XY,1,rank_XY(),sendtag+7);
req2[7] = Comm.Irecv(&recvCount_xy,1,rank_xy(),recvtag+7);
req1[8] = Comm.Isend(&sendCount_Xy,1,rank_Xy(),sendtag+8);
req2[8] = Comm.Irecv(&recvCount_xY,1,rank_xY(),recvtag+8);
req1[9] = Comm.Isend(&sendCount_xY,1,rank_xY(),sendtag+9);
req2[9] = Comm.Irecv(&recvCount_Xy,1,rank_Xy(),recvtag+9);
req1[10] = Comm.Isend(&sendCount_xz,1,rank_xz(),sendtag+10);
req2[10] = Comm.Irecv(&recvCount_XZ,1,rank_XZ(),recvtag+10);
req1[11] = Comm.Isend(&sendCount_XZ,1,rank_XZ(),sendtag+11);
req2[11] = Comm.Irecv(&recvCount_xz,1,rank_xz(),recvtag+11);
req1[12] = Comm.Isend(&sendCount_Xz,1,rank_Xz(),sendtag+12);
req2[12] = Comm.Irecv(&recvCount_xZ,1,rank_xZ(),recvtag+12);
req1[13] = Comm.Isend(&sendCount_xZ,1,rank_xZ(),sendtag+13);
req2[13] = Comm.Irecv(&recvCount_Xz,1,rank_Xz(),recvtag+13);
req1[14] = Comm.Isend(&sendCount_yz,1,rank_yz(),sendtag+14);
req2[14] = Comm.Irecv(&recvCount_YZ,1,rank_YZ(),recvtag+14);
req1[15] = Comm.Isend(&sendCount_YZ,1,rank_YZ(),sendtag+15);
req2[15] = Comm.Irecv(&recvCount_yz,1,rank_yz(),recvtag+15);
req1[16] = Comm.Isend(&sendCount_Yz,1,rank_Yz(),sendtag+16);
req2[16] = Comm.Irecv(&recvCount_yZ,1,rank_yZ(),recvtag+16);
req1[17] = Comm.Isend(&sendCount_yZ,1,rank_yZ(),sendtag+17);
req2[17] = Comm.Irecv(&recvCount_Yz,1,rank_Yz(),recvtag+17);
Comm.waitAll(18,req1);
Comm.waitAll(18,req2);
Comm.barrier();
// allocate recv lists
recvList_x.resize( recvCount_x, 0 );
recvList_y.resize( recvCount_y, 0 );
recvList_z.resize( recvCount_z, 0 );
recvList_X.resize( recvCount_X, 0 );
recvList_Y.resize( recvCount_Y, 0 );
recvList_Z.resize( recvCount_Z, 0 );
recvList_xy.resize( recvCount_xy, 0 );
recvList_yz.resize( recvCount_yz, 0 );
recvList_xz.resize( recvCount_xz, 0 );
recvList_Xy.resize( recvCount_Xy, 0 );
recvList_Yz.resize( recvCount_Yz, 0 );
recvList_xZ.resize( recvCount_xZ, 0 );
recvList_xY.resize( recvCount_xY, 0 );
recvList_yZ.resize( recvCount_yZ, 0 );
recvList_Xz.resize( recvCount_Xz, 0 );
recvList_XY.resize( recvCount_XY, 0 );
recvList_YZ.resize( recvCount_YZ, 0 );
recvList_XZ.resize( recvCount_XZ, 0 );
//......................................................................................
req1[0] = Comm.Isend(sendList_x.data(),sendCount_x,rank_x(),sendtag);
req2[0] = Comm.Irecv(recvList_X.data(),recvCount_X,rank_X(),recvtag);
req1[1] = Comm.Isend(sendList_X.data(),sendCount_X,rank_X(),sendtag);
req2[1] = Comm.Irecv(recvList_x.data(),recvCount_x,rank_x(),recvtag);
req1[2] = Comm.Isend(sendList_y.data(),sendCount_y,rank_y(),sendtag);
req2[2] = Comm.Irecv(recvList_Y.data(),recvCount_Y,rank_Y(),recvtag);
req1[3] = Comm.Isend(sendList_Y.data(),sendCount_Y,rank_Y(),sendtag);
req2[3] = Comm.Irecv(recvList_y.data(),recvCount_y,rank_y(),recvtag);
req1[4] = Comm.Isend(sendList_z.data(),sendCount_z,rank_z(),sendtag);
req2[4] = Comm.Irecv(recvList_Z.data(),recvCount_Z,rank_Z(),recvtag);
req1[5] = Comm.Isend(sendList_Z.data(),sendCount_Z,rank_Z(),sendtag);
req2[5] = Comm.Irecv(recvList_z.data(),recvCount_z,rank_z(),recvtag);
req1[6] = Comm.Isend(sendList_xy.data(),sendCount_xy,rank_xy(),sendtag);
req2[6] = Comm.Irecv(recvList_XY.data(),recvCount_XY,rank_XY(),recvtag);
req1[7] = Comm.Isend(sendList_XY.data(),sendCount_XY,rank_XY(),sendtag);
req2[7] = Comm.Irecv(recvList_xy.data(),recvCount_xy,rank_xy(),recvtag);
req1[8] = Comm.Isend(sendList_Xy.data(),sendCount_Xy,rank_Xy(),sendtag);
req2[8] = Comm.Irecv(recvList_xY.data(),recvCount_xY,rank_xY(),recvtag);
req1[9] = Comm.Isend(sendList_xY.data(),sendCount_xY,rank_xY(),sendtag);
req2[9] = Comm.Irecv(recvList_Xy.data(),recvCount_Xy,rank_Xy(),recvtag);
req1[10] = Comm.Isend(sendList_xz.data(),sendCount_xz,rank_xz(),sendtag);
req2[10] = Comm.Irecv(recvList_XZ.data(),recvCount_XZ,rank_XZ(),recvtag);
req1[11] = Comm.Isend(sendList_XZ.data(),sendCount_XZ,rank_XZ(),sendtag);
req2[11] = Comm.Irecv(recvList_xz.data(),recvCount_xz,rank_xz(),recvtag);
req1[12] = Comm.Isend(sendList_Xz.data(),sendCount_Xz,rank_Xz(),sendtag);
req2[12] = Comm.Irecv(recvList_xZ.data(),recvCount_xZ,rank_xZ(),recvtag);
req1[13] = Comm.Isend(sendList_xZ.data(),sendCount_xZ,rank_xZ(),sendtag);
req2[13] = Comm.Irecv(recvList_Xz.data(),recvCount_Xz,rank_Xz(),recvtag);
req1[14] = Comm.Isend(sendList_yz.data(),sendCount_yz,rank_yz(),sendtag);
req2[14] = Comm.Irecv(recvList_YZ.data(),recvCount_YZ,rank_YZ(),recvtag);
req1[15] = Comm.Isend(sendList_YZ.data(),sendCount_YZ,rank_YZ(),sendtag);
req2[15] = Comm.Irecv(recvList_yz.data(),recvCount_yz,rank_yz(),recvtag);
req1[16] = Comm.Isend(sendList_Yz.data(),sendCount_Yz,rank_Yz(),sendtag);
req2[16] = Comm.Irecv(recvList_yZ.data(),recvCount_yZ,rank_yZ(),recvtag);
req1[17] = Comm.Isend(sendList_yZ.data(),sendCount_yZ,rank_yZ(),sendtag);
req2[17] = Comm.Irecv(recvList_Yz.data(),recvCount_Yz,rank_Yz(),recvtag);
Comm.waitAll(18,req1);
Comm.waitAll(18,req2);
//......................................................................................
for (int idx=0; idx<recvCount_x; idx++) recvList_x[idx] -= (Nx-2);
for (int idx=0; idx<recvCount_X; idx++) recvList_X[idx] += (Nx-2);
for (int idx=0; idx<recvCount_y; idx++) recvList_y[idx] -= (Ny-2)*Nx;
for (int idx=0; idx<recvCount_Y; idx++) recvList_Y[idx] += (Ny-2)*Nx;
for (int idx=0; idx<recvCount_z; idx++) recvList_z[idx] -= (Nz-2)*Nx*Ny;
for (int idx=0; idx<recvCount_Z; idx++) recvList_Z[idx] += (Nz-2)*Nx*Ny;
for (int idx=0; idx<recvCount_xy; idx++) recvList_xy[idx] -= (Nx-2)+(Ny-2)*Nx;
for (int idx=0; idx<recvCount_XY; idx++) recvList_XY[idx] += (Nx-2)+(Ny-2)*Nx;
for (int idx=0; idx<recvCount_xY; idx++) recvList_xY[idx] -= (Nx-2)-(Ny-2)*Nx;
for (int idx=0; idx<recvCount_Xy; idx++) recvList_Xy[idx] += (Nx-2)-(Ny-2)*Nx;
for (int idx=0; idx<recvCount_xz; idx++) recvList_xz[idx] -= (Nx-2)+(Nz-2)*Nx*Ny;
for (int idx=0; idx<recvCount_XZ; idx++) recvList_XZ[idx] += (Nx-2)+(Nz-2)*Nx*Ny;
for (int idx=0; idx<recvCount_xZ; idx++) recvList_xZ[idx] -= (Nx-2)-(Nz-2)*Nx*Ny;
for (int idx=0; idx<recvCount_Xz; idx++) recvList_Xz[idx] += (Nx-2)-(Nz-2)*Nx*Ny;
for (int idx=0; idx<recvCount_yz; idx++) recvList_yz[idx] -= (Ny-2)*Nx + (Nz-2)*Nx*Ny;
for (int idx=0; idx<recvCount_YZ; idx++) recvList_YZ[idx] += (Ny-2)*Nx + (Nz-2)*Nx*Ny;
for (int idx=0; idx<recvCount_yZ; idx++) recvList_yZ[idx] -= (Ny-2)*Nx - (Nz-2)*Nx*Ny;
for (int idx=0; idx<recvCount_Yz; idx++) recvList_Yz[idx] += (Ny-2)*Nx - (Nz-2)*Nx*Ny;
//......................................................................................
//......................................................................................
}
void Domain::ReadIDs(){
// Read the IDs from input file
int nprocs=nprocx()*nprocy()*nprocz();
size_t readID;
char LocalRankString[8];
char LocalRankFilename[40];
//.......................................................................
if (rank() == 0) printf("Read input media... \n");
//.......................................................................
sprintf(LocalRankString,"%05d",rank());
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// .......... READ THE INPUT FILE .......................................
if (rank()==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
sprintf(LocalRankFilename,"ID.%05i",rank());
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (!IDFILE) ERROR("Domain::ReadIDs -- Error opening file: ID.xxxxx");
readID=fread(id.data(),1,N,IDFILE);
if (readID != size_t(N)) printf("Domain::ReadIDs -- Error reading ID (rank=%i) \n",rank());
fclose(IDFILE);
// Compute the porosity
double sum;
double sum_local=0.0;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
//.........................................................
// If external boundary conditions are applied remove solid
if (BoundaryCondition > 0 && kproc() == 0){
if (inlet_layers_z < 4) inlet_layers_z=4;
for (int k=0; k<inlet_layers_z; k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
id[n] = 1;
}
}
}
}
if (BoundaryCondition > 0 && kproc() == nprocz()-1){
if (outlet_layers_z < 4) outlet_layers_z=4;
for (int k=Nz-outlet_layers_z; 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;
id[n] = 2;
}
}
}
}
for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-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;
if (id[n] > 0){
sum_local+=1.0;
}
}
}
}
sum = Comm.sumReduce(sum_local);
porosity = sum*iVol_global;
if (rank()==0) printf("Media porosity = %f \n",porosity);
//.........................................................
}
int Domain::PoreCount(){
/*
* count the number of nodes occupied by mobile phases
*/
int Npore=0; // number of local pore nodes
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;
if (id[n] > 0){
Npore++;
}
}
}
}
return Npore;
}
void Domain::CommunicateMeshHalo(DoubleArray &Mesh)
{
int sendtag, recvtag;
sendtag = recvtag = 7;
double *MeshData = Mesh.data();
// send buffers
auto sendData_x = new double [sendCount("x")];
auto sendData_y = new double [sendCount("y")];
auto sendData_z = new double [sendCount("z")];
auto sendData_X = new double [sendCount("X")];
auto sendData_Y = new double [sendCount("Y")];
auto sendData_Z = new double [sendCount("Z")];
auto sendData_xy = new double [sendCount("xy")];
auto sendData_yz = new double [sendCount("yz")];
auto sendData_xz = new double [sendCount("xz")];
auto sendData_Xy = new double [sendCount("Xy")];
auto sendData_Yz = new double [sendCount("Yz")];
auto sendData_xZ = new double [sendCount("xZ")];
auto sendData_xY = new double [sendCount("xY")];
auto sendData_yZ = new double [sendCount("yZ")];
auto sendData_Xz = new double [sendCount("Xz")];
auto sendData_XY = new double [sendCount("XY")];
auto sendData_YZ = new double [sendCount("YZ")];
auto sendData_XZ = new double [sendCount("XZ")];
// recv buffers
auto recvData_x = new double [recvCount("x")];
auto recvData_y = new double [recvCount("y")];
auto recvData_z = new double [recvCount("z")];
auto recvData_X = new double [recvCount("X")];
auto recvData_Y = new double [recvCount("Y")];
auto recvData_Z = new double [recvCount("Z")];
auto recvData_xy = new double [recvCount("xy")];
auto recvData_yz = new double [recvCount("yz")];
auto recvData_xz = new double [recvCount("xz")];
auto recvData_Xy = new double [recvCount("Xy")];
auto recvData_xZ = new double [recvCount("xZ")];
auto recvData_xY = new double [recvCount("xY")];
auto recvData_yZ = new double [recvCount("yZ")];
auto recvData_Yz = new double [recvCount("Yz")];
auto recvData_Xz = new double [recvCount("Xz")];
auto recvData_XY = new double [recvCount("XY")];
auto recvData_YZ = new double [recvCount("YZ")];
auto recvData_XZ = new double [recvCount("XZ")];
// Pack data
PackMeshData(sendList("x"), sendCount("x"), sendData_x, MeshData);
PackMeshData(sendList("X"), sendCount("X"), sendData_X, MeshData);
PackMeshData(sendList("y"), sendCount("y"), sendData_y, MeshData);
PackMeshData(sendList("Y"), sendCount("Y"), sendData_Y, MeshData);
PackMeshData(sendList("z"), sendCount("z"), sendData_z, MeshData);
PackMeshData(sendList("Z"), sendCount("Z"), sendData_Z, MeshData);
PackMeshData(sendList("xy"), sendCount("xy"), sendData_xy, MeshData);
PackMeshData(sendList("Xy"), sendCount("Xy"), sendData_Xy, MeshData);
PackMeshData(sendList("xY"), sendCount("xY"), sendData_xY, MeshData);
PackMeshData(sendList("XY"), sendCount("XY"), sendData_XY, MeshData);
PackMeshData(sendList("xz"), sendCount("xz"), sendData_xz, MeshData);
PackMeshData(sendList("Xz"), sendCount("Xz"), sendData_Xz, MeshData);
PackMeshData(sendList("xZ"), sendCount("xZ"), sendData_xZ, MeshData);
PackMeshData(sendList("XZ"), sendCount("XZ"), sendData_XZ, MeshData);
PackMeshData(sendList("yz"), sendCount("yz"), sendData_yz, MeshData);
PackMeshData(sendList("Yz"), sendCount("Yz"), sendData_Yz, MeshData);
PackMeshData(sendList("yZ"), sendCount("yZ"), sendData_yZ, MeshData);
PackMeshData(sendList("YZ"), sendCount("YZ"), sendData_YZ, MeshData);
// send/recv
Comm.sendrecv(sendData_x,sendCount("x"),rank_x(),sendtag,recvData_X,recvCount("X"),rank_X(),recvtag);
Comm.sendrecv(sendData_X,sendCount("X"),rank_X(),sendtag,recvData_x,recvCount("x"),rank_x(),recvtag);
Comm.sendrecv(sendData_y,sendCount("y"),rank_y(),sendtag,recvData_Y,recvCount("Y"),rank_Y(),recvtag);
Comm.sendrecv(sendData_Y,sendCount("Y"),rank_Y(),sendtag,recvData_y,recvCount("y"),rank_y(),recvtag);
Comm.sendrecv(sendData_z,sendCount("z"),rank_z(),sendtag,recvData_Z,recvCount("Z"),rank_Z(),recvtag);
Comm.sendrecv(sendData_Z,sendCount("Z"),rank_Z(),sendtag,recvData_z,recvCount("z"),rank_z(),recvtag);
Comm.sendrecv(sendData_xy,sendCount("xy"),rank_xy(),sendtag,recvData_XY,recvCount("XY"),rank_XY(),recvtag);
Comm.sendrecv(sendData_XY,sendCount("XY"),rank_XY(),sendtag,recvData_xy,recvCount("xy"),rank_xy(),recvtag);
Comm.sendrecv(sendData_Xy,sendCount("Xy"),rank_Xy(),sendtag,recvData_xY,recvCount("xY"),rank_xY(),recvtag);
Comm.sendrecv(sendData_xY,sendCount("xY"),rank_xY(),sendtag,recvData_Xy,recvCount("Xy"),rank_Xy(),recvtag);
Comm.sendrecv(sendData_xz,sendCount("xz"),rank_xz(),sendtag,recvData_XZ,recvCount("XZ"),rank_XZ(),recvtag);
Comm.sendrecv(sendData_XZ,sendCount("XZ"),rank_XZ(),sendtag,recvData_xz,recvCount("xz"),rank_xz(),recvtag);
Comm.sendrecv(sendData_Xz,sendCount("Xz"),rank_Xz(),sendtag,recvData_xZ,recvCount("xZ"),rank_xZ(),recvtag);
Comm.sendrecv(sendData_xZ,sendCount("xZ"),rank_xZ(),sendtag,recvData_Xz,recvCount("Xz"),rank_Xz(),recvtag);
Comm.sendrecv(sendData_yz,sendCount("yz"),rank_yz(),sendtag,recvData_YZ,recvCount("YZ"),rank_YZ(),recvtag);
Comm.sendrecv(sendData_YZ,sendCount("YZ"),rank_YZ(),sendtag,recvData_yz,recvCount("yz"),rank_yz(),recvtag);
Comm.sendrecv(sendData_Yz,sendCount("Yz"),rank_Yz(),sendtag,recvData_yZ,recvCount("yZ"),rank_yZ(),recvtag);
Comm.sendrecv(sendData_yZ,sendCount("yZ"),rank_yZ(),sendtag,recvData_Yz,recvCount("Yz"),rank_Yz(),recvtag);
// unpack data
UnpackMeshData(recvList("x"), recvCount("x") ,recvData_x, MeshData);
UnpackMeshData(recvList("X"), recvCount("X") ,recvData_X, MeshData);
UnpackMeshData(recvList("y"), recvCount("y") ,recvData_y, MeshData);
UnpackMeshData(recvList("Y"), recvCount("Y") ,recvData_Y, MeshData);
UnpackMeshData(recvList("z"), recvCount("z") ,recvData_z, MeshData);
UnpackMeshData(recvList("Z"), recvCount("Z") ,recvData_Z, MeshData);
UnpackMeshData(recvList("xy"), recvCount("xy") ,recvData_xy, MeshData);
UnpackMeshData(recvList("Xy"), recvCount("Xy") ,recvData_Xy, MeshData);
UnpackMeshData(recvList("xY"), recvCount("xY") ,recvData_xY, MeshData);
UnpackMeshData(recvList("XY"), recvCount("XY") ,recvData_XY, MeshData);
UnpackMeshData(recvList("xz"), recvCount("xz") ,recvData_xz, MeshData);
UnpackMeshData(recvList("Xz"), recvCount("Xz") ,recvData_Xz, MeshData);
UnpackMeshData(recvList("xZ"), recvCount("xZ") ,recvData_xZ, MeshData);
UnpackMeshData(recvList("XZ"), recvCount("XZ") ,recvData_XZ, MeshData);
UnpackMeshData(recvList("yz"), recvCount("yz") ,recvData_yz, MeshData);
UnpackMeshData(recvList("Yz"), recvCount("Yz") ,recvData_Yz, MeshData);
UnpackMeshData(recvList("yZ"), recvCount("yZ") ,recvData_yZ, MeshData);
UnpackMeshData(recvList("YZ"), recvCount("YZ") ,recvData_YZ, MeshData);
// 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;
}
// Ideally stuff below here should be moved somewhere else -- doesn't really belong here
void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, size_t Np)
{
double value;
ofstream File(FILENAME,ios::binary);
for (size_t n=0; n<Np; n++){
// Write the two density values
value = cDen[n];
File.write((char*) &value, sizeof(value));
value = cDen[Np+n];
File.write((char*) &value, sizeof(value));
// Write the even distributions
for (size_t q=0; q<19; q++){
value = cfq[q*Np+n];
File.write((char*) &value, sizeof(value));
}
}
File.close();
}
void ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, size_t Np)
{
double value=0;
ifstream File(FILENAME,ios::binary);
for (size_t n=0; n<Np; n++){
File.read((char*) &value, sizeof(value));
cPhi[n] = value;
// Read the distributions
for (size_t q=0; q<19; q++){
File.read((char*) &value, sizeof(value));
cfq[q*Np+n] = value;
}
}
File.close();
}
void ReadBinaryFile(char *FILENAME, double *Data, size_t N)
{
double value;
ifstream File(FILENAME,ios::binary);
if (File.good()){
for (size_t n=0; n<N; n++){
// Write the two density values
File.read((char*) &value, sizeof(value));
Data[n] = value;
}
}
else {
for (size_t n=0; n<N; n++) Data[n] = 1.2e-34;
}
File.close();
}
void Domain::ReadFromFile(const std::string& Filename,const std::string& Datatype, double *UserData)
{
//........................................................................................
// Reading the user-defined input file
// NOTE: so far it only supports BC=0 (periodic) and BC=5 (mixed reflection)
// because if checkerboard or inlet/outlet buffer layers are added, the
// value of the void space is undefined.
// NOTE: if BC=5 is used, where the inlet and outlet layers of the domain are modified,
// user needs to modify the input file accordingly before LBPM simulator read
// the input file.
//........................................................................................
int RANK = rank();
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz;
int64_t global_Nx,global_Ny,global_Nz;
int64_t i,j,k;
//TODO These offset we may still need them
int64_t xStart,yStart,zStart;
xStart=yStart=zStart=0;
// Read domain parameters
// TODO currently the size of the data is still read from Domain{};
// but user may have a user-specified size
auto size = database->getVector<int>( "n" );
auto SIZE = database->getVector<int>( "N" );
auto nproc = database->getVector<int>( "nproc" );
//TODO currently the funcationality "offset" is disabled as the user-defined input data may have a different size from that of the input domain
if (database->keyExists( "offset" )){
auto offset = database->getVector<int>( "offset" );
xStart = offset[0];
yStart = offset[1];
zStart = offset[2];
}
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;
double *SegData = NULL;
if (RANK==0){
printf("User-defined input file: %s (data type: %s)\n",Filename.c_str(),Datatype.c_str());
printf("NOTE: currently only BC=0 or 5 supports user-defined input file!\n");
// Rank=0 reads the entire segmented data and distributes to worker processes
printf("Dimensions of the user-defined input file: %ld x %ld x %ld \n",global_Nx,global_Ny,global_Nz);
int64_t SIZE = global_Nx*global_Ny*global_Nz;
if (Datatype == "double"){
printf("Reading input data as double precision floating number\n");
SegData = new double[SIZE];
FILE *SEGDAT = fopen(Filename.c_str(),"rb");
if (SEGDAT==NULL) ERROR("Domain.cpp: Error reading user-defined file!\n");
size_t ReadSeg;
ReadSeg=fread(SegData,8,SIZE,SEGDAT);
if (ReadSeg != size_t(SIZE)) printf("Domain.cpp: Error reading file: %s\n",Filename.c_str());
fclose(SEGDAT);
}
else{
ERROR("Error: User-defined input file only supports double-precision floating number!\n");
}
printf("Read file successfully from %s \n",Filename.c_str());
}
// 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;
int64_t z_transition_size = 0;
//char LocalRankFilename[1000];//just for debug
double *loc_id;
loc_id = new double [(nx+2)*(ny+2)*(nz+2)];
// Set up the sub-domains
if (RANK==0){
printf("Decomposing user-defined input file\n");
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];
}
}
}
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;
UserData[nlocal] = loc_id[nlocal];
}
}
}
}
else{
//printf("Sending data to process %i \n", rnk);
Comm.send(loc_id,N,rnk,15);
}
}
}
}
}
else{
// Recieve the subdomain from rank = 0
//printf("Ready to recieve data %i at process %i \n", N,rank);
Comm.recv(id.data(),N,0,15);
}
Comm.barrier();
}
void Domain::AggregateLabels( const std::string& filename, DoubleArray &UserData ){
int nx = Nx;
int ny = Ny;
int nz = Nz;
int npx = nprocx();
int npy = nprocy();
int npz = nprocz();
int ipx = iproc();
int ipy = jproc();
int ipz = kproc();
int nprocs = nprocx()*nprocy()*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);
unsigned long int full_size = long(full_nx)*long(full_ny)*long(full_nz);
double *LocalID;
LocalID = new double [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++){
double local_id_val = UserData(i,j,k);
LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val;
}
}
}
Comm.barrier();
// populate the FullID
if (rank() == 0){
double *FullID;
FullID = new double [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;
unsigned long int n_full = z*long(full_nx)*long(full_ny) + y*long(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_DOUBLE,rnk,tag,Comm,MPI_STATUS_IGNORE);
Comm.recv(LocalID,local_size,rnk,tag);
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;
unsigned long int n_full = z*long(full_nx)*long(full_ny) + y*long(full_nx) + x;
FullID[n_full] = LocalID[n_local];
}
}
}
}
// write the output
FILE *OUTFILE = fopen(filename.c_str(),"wb");
fwrite(FullID,8,full_size,OUTFILE);
fclose(OUTFILE);
}
else{
// send LocalID to rank=0
int tag = 15+ rank();
int dstrank = 0;
//MPI_Send(LocalID,local_size,MPI_DOUBLE,dstrank,tag,Comm);
Comm.send(LocalID,local_size,dstrank,tag);
}
Comm.barrier();
}
Reference in New Issue View Git Blame Copy Permalink