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

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

Browse Source
This commit is contained in:
James E McClure 2018-06-11 14:57:47 -04:00
commit 8bd361ae03
14 changed files with 930 additions and 435 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
IO/netcdf.cpp
View File

@ -317,7 +317,7 @@ Array<std::string> getVar<std::string>( int fid, const std::string& var )
{
PROFILE_START("getVar<std::string>");
Array<char> tmp = getVar<char>( fid, var );
std::vector<size_t> dim = tmp.size();
std::vector<size_t> dim = {tmp.size(0), tmp.size(1), tmp.size(2) };
if ( dim.size() == 1 )
dim[0] = 1;
else
@ -451,7 +451,7 @@ void write( int fid, const std::string& var, const std::vector<int>& dimids,
CHECK_NC_ERR( err );
// parallel write: each process writes its subarray to the file
auto x = data.reverseDim();
std::vector<size_t> count = data.size();
std::vector<size_t> count = { data.size(0), data.size(1), data.size(2) };
std::vector<size_t> start = { info.ix*data.size(0), info.jy*data.size(1), info.kz*data.size(2) };
nc_put_vara( fid, varid, start.data(), count.data(), x.data() );
}

227
analysis/Minkowski.cpp Normal file
View File

@ -0,0 +1,227 @@
#include "analysis/Minkowski.h"
#include "analysis/pmmc.h"
#include "common/Domain.h"
#include "common/Communication.h"
#include "analysis/analysis.h"
#include "shared_ptr.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "IO/MeshDatabase.h"
#include "IO/Reader.h"
#include "IO/Writer.h"
#define PI 3.14159265359
// Constructor
Minkowski::Minkowski(std::shared_ptr <Domain> dm):
n_obj_pts(0), n_obj_tris(0), kstart(0), kfinish(0), isovalue(0), Volume(0),
LOGFILE(NULL), Dm(dm), vol_n(0), vol_n_global(0)
{
Nx=dm->Nx; Ny=dm->Ny; Nz=dm->Nz;
Volume=double((Nx-2)*(Ny-2)*(Nz-2))*double(Dm->nprocx()*Dm->nprocy()*Dm->nprocz());
TempID = new char[Nx*Ny*Nz];
// Global arrays
PhaseID.resize(Nx,Ny,Nz); PhaseID.fill(0);
SDn.resize(Nx,Ny,Nz); SDn.fill(0);
MeanCurvature.resize(Nx,Ny,Nz); MeanCurvature.fill(0);
GaussCurvature.resize(Nx,Ny,Nz); GaussCurvature.fill(0);
SDn_x.resize(Nx,Ny,Nz); SDn_x.fill(0); // Gradient of the signed distance
SDn_y.resize(Nx,Ny,Nz); SDn_y.fill(0);
SDn_z.resize(Nx,Ny,Nz); SDn_z.fill(0);
//.........................................
// Allocate cube storage space
CubeValues.resize(2,2,2);
obj_tris.resize(3,20);
obj_pts=DTMutableList<Point>(20);
tmp=DTMutableList<Point>(20);
//.........................................
Values.resize(20);
//DistanceValues.resize(20);
NormalVector.resize(60);
if (Dm->rank()==0){
LOGFILE = fopen("minkowski.csv","a+");
if (fseek(LOGFILE,0,SEEK_SET) == fseek(LOGFILE,0,SEEK_CUR))
{
// If LOGFILE is empty, write a short header to list the averages
//fprintf(LOGFILE,"--------------------------------------------------------------------------------------\n");
fprintf(LOGFILE,"Vn An Jn Xn\n"); //miknowski measures,
}
}
}
// Destructor
Minkowski::~Minkowski()
{
if ( LOGFILE!=NULL ) { fclose(LOGFILE); }
}
void Minkowski::Initialize()
{
isovalue=0.0;
vol_n = euler = Jn = An = Kn = 0.0;
}
void Minkowski::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);
//...........................................................................
// Gradient of the phase indicator field
//...........................................................................
Dm->CommunicateMeshHalo(SDn_x);
//...........................................................................
Dm->CommunicateMeshHalo(SDn_y);
//...........................................................................
Dm->CommunicateMeshHalo(SDn_z);
//...........................................................................
//...........................................................................
// Compute the mesh curvature of the phase indicator field
pmmc_MeshCurvature(SDn, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
//...........................................................................
//...........................................................................
Dm->CommunicateMeshHalo(MeanCurvature);
//...........................................................................
Dm->CommunicateMeshHalo(GaussCurvature);
//...........................................................................
// Initializing the blob ID
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 (Dm->id[n] == 0){
// Solid phase
PhaseID(i,j,k) = 0;
}
else {
// non-wetting phase
PhaseID(i,j,k) = 1;
}
}
}
}
}
void Minkowski::ComputeLocal()
{
int i,j,k,n,kmin,kmax;
int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};
// If external boundary conditions are set, do not average over the inlet
kmin=1; kmax=Nz-1;
if (Dm->BoundaryCondition > 0 && Dm->kproc() == 0) kmin=4;
if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
vol_n = euler = Jn = An = Kn = 0.0;
for (k=kmin; k<kmax; k++){
for (j=1; j<Ny-1; j++){
for (i=1; i<Nx-1; i++){
//...........................................................................
n_obj_pts=0;
n_obj_tris=0;
//...........................................................................
// Compute volume averages
for (int p=0;p<8;p++){
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
if ( Dm->id[n] != 0 ){
// 1-D index for this cube corner
if ( SDn(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < 0 ){
vol_n += 0.125;
}
}
}
n_obj_pts=n_obj_tris=0;
// Compute the non-wetting phase surface and associated area
An += geomavg_MarchingCubes(SDn,isovalue,i,j,k,obj_pts,n_obj_pts,obj_tris,n_obj_tris);
Jn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,obj_pts,obj_tris,Values,
i,j,k,n_obj_pts,n_obj_tris);
// Compute Euler characteristic from integral of gaussian curvature
Kn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,obj_pts,obj_tris,Values,
i,j,k,n_obj_pts,n_obj_tris);
euler += geomavg_EulerCharacteristic(obj_pts,obj_tris,n_obj_pts,n_obj_tris,i,j,k);
}
}
}
}
/*
void Minkowski::AssignComponentLabels()
{
int LabelNWP=1;
int LabelWP=2;
// NOTE: labeling the wetting phase components is tricky! One sandstone media had over 800,000 components
// NumberComponents_WP = ComputeGlobalPhaseComponent(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->rank_info,PhaseID,LabelWP,Label_WP);
// treat all wetting phase is connected
NumberComponents_WP=1;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
Label_WP(i,j,k) = 0;
//if (SDs(i,j,k) > 0.0) PhaseID(i,j,k) = 0;
//else if (Phase(i,j,k) > 0.0) PhaseID(i,j,k) = LabelNWP;
//else PhaseID(i,j,k) = LabelWP;
}
}
}
// Fewer non-wetting phase features are present
//NumberComponents_NWP = ComputeGlobalPhaseComponent(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->rank_info,PhaseID,LabelNWP,Label_NWP);
NumberComponents_NWP = ComputeGlobalBlobIDs(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->rank_info,SDs,SDn,solid_isovalue,isovalue,Label_NWP,Dm->Comm);
}
*/
void Minkowski::Reduce()
{
int i;
double iVol_global=1.0/Volume;
//...........................................................................
MPI_Barrier(Dm->Comm);
// Phase averages
MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&euler,&euler_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&An,&An_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Allreduce(&Jn,&Jn_global,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
MPI_Barrier(Dm->Comm);
// normalize to per unit volume basis
vol_n_global /=Volume;
An_global /=Volume;
Jn_global /=Volume;
euler_global /=Volume;
euler_global /= (2*PI);
}
void Minkowski::NonDimensionalize(double D)
{
An_global *= D;
Jn_global *= D*D;
euler_global *= D*D*D;
}
void Minkowski::PrintAll()
{
if (Dm->rank()==0){
fprintf(LOGFILE,"%.5g %.5g %.5g %.5g\n",vol_n_global, An_global, Jn_global, euler_global); // minkowski measures
fflush(LOGFILE);
}
}

87
analysis/Minkowski.h Normal file
View File

@ -0,0 +1,87 @@
// Header file for two-phase averaging class
#ifndef Minkowski_INC
#define Minkowski_INC
#include <vector>
#include "analysis/pmmc.h"
#include "common/Domain.h"
#include "common/Communication.h"
#include "analysis/analysis.h"
#include "shared_ptr.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "IO/MeshDatabase.h"
#include "IO/Reader.h"
#include "IO/Writer.h"
class Minkowski{
//...........................................................................
int n_obj_pts;
int n_obj_tris;
//...........................................................................
int nc;
int kstart,kfinish;
double isovalue;
double Volume;
// initialize lists for vertices for surfaces, common line
DTMutableList<Point> obj_pts;
DTMutableList<Point> tmp;
// initialize triangle lists for surfaces
IntArray obj_tris;
// Temporary storage arrays
DoubleArray CubeValues;
DoubleArray Values;
DoubleArray NormalVector;
DoubleArray RecvBuffer;
char *TempID;
// CSV / text file where time history of averages is saved
FILE *LOGFILE;
public:
//...........................................................................
std::shared_ptr <Domain> Dm;
//...........................................................................
// Averaging variables
//...........................................................................
// local averages (to each MPI process)
double vol_n; // volumes the exclude the interfacial region
// Global averages (all processes)
double vol_n_global; // volumes the exclude the interfacial region
double euler,Kn,Jn,An;
double euler_global,Kn_global,Jn_global,An_global;
//...........................................................................
int Nx,Ny,Nz;
IntArray PhaseID; // Phase ID array (solid=0, non-wetting=1, wetting=2)
DoubleArray SDn;
DoubleArray MeanCurvature;
DoubleArray GaussCurvature;
DoubleArray SDn_x; // Gradient of the signed distance
DoubleArray SDn_y;
DoubleArray SDn_z;
//...........................................................................
Minkowski(std::shared_ptr <Domain> Dm);
~Minkowski();
void Initialize();
void UpdateMeshValues();
void ComputeLocal();
void Reduce();
void NonDimensionalize(double D);
void PrintAll();
int GetCubeLabel(int i, int j, int k, IntArray &BlobLabel);
void SortBlobs();
};
#endif

20
analysis/uCT.cpp
View File

@ -141,7 +141,8 @@ void removeDisconnected( Array<char>& ID, const Domain& Dm )
// Solve a level (without any coarse level information)
void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
Array<float>& Dist, Array<float>& MultiScaleSmooth, Array<float>& NonLocalMean,
fillHalo<float>& fillFloat, const Domain& Dm, int nprocx )
fillHalo<float>& fillFloat, const Domain& Dm, int nprocx,
float threshold, float lamda, float sigsq, int depth)
{
PROFILE_SCOPED(timer,"solve");
// Compute the median filter on the sparse array
@ -153,8 +154,8 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
fillFloat.fill(Dist);
smooth( VOL, Dist, 2.0, MultiScaleSmooth, fillFloat );
// Compute non-local mean
int depth = 5;
float sigsq=0.1;
// int depth = 5;
// float sigsq=0.1;
int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq);
fillFloat.fill(NonLocalMean);
}
@ -164,7 +165,8 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
void refine( const Array<float>& Dist_coarse,
const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
Array<float>& Dist, Array<float>& MultiScaleSmooth, Array<float>& NonLocalMean,
fillHalo<float>& fillFloat, const Domain& Dm, int nprocx, int level )
fillHalo<float>& fillFloat, const Domain& Dm, int nprocx, int level,
float threshold, float lamda, float sigsq, int depth)
{
PROFILE_SCOPED(timer,"refine");
int ratio[3] = { int(Dist.size(0)/Dist_coarse.size(0)),
@ -175,7 +177,7 @@ void refine( const Array<float>& Dist_coarse,
// Compute the median filter on the array and segment
Med3D( VOL, Mean );
fillFloat.fill( Mean );
segment( Mean, ID, 0.01 );
segment( Mean, ID, threshold );
// If the ID has the wrong distance, set the distance to 0 and run a simple filter to set neighbors to 0
for (size_t i=0; i<ID.length(); i++) {
char id = Dist(i)>0 ? 1:0;
@ -193,11 +195,11 @@ void refine( const Array<float>& Dist_coarse,
Dist = imfilter::imfilter_separable<float>( Dist, {1,1,1}, filter_set, BC );
fillFloat.fill( Dist );
// Smooth the volume data
float lambda = 2*sqrt(double(ratio[0]*ratio[0]+ratio[1]*ratio[1]+ratio[2]*ratio[2]));
smooth( VOL, Dist, lambda, MultiScaleSmooth, fillFloat );
float h = 2*lamda*sqrt(double(ratio[0]*ratio[0]+ratio[1]*ratio[1]+ratio[2]*ratio[2]));
smooth( VOL, Dist, h, MultiScaleSmooth, fillFloat );
// Compute non-local mean
int depth = 3;
float sigsq = 0.1;
// int depth = 3;
// float sigsq = 0.1;
int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq);
fillFloat.fill(NonLocalMean);
segment( NonLocalMean, ID, 0.001 );

6
analysis/uCT.h
View File

@ -31,14 +31,16 @@ void removeDisconnected( Array<char>& ID, const Domain& Dm );
// Solve a level (without any coarse level information)
void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
Array<float>& Dist, Array<float>& MultiScaleSmooth, Array<float>& NonLocalMean,
fillHalo<float>& fillFloat, const Domain& Dm, int nprocx );
fillHalo<float>& fillFloat, const Domain& Dm, int nprocx,
float threshold, float lamda, float sigsq, int depth);
// Refine a solution from a coarse grid to a fine grid
void refine( const Array<float>& Dist_coarse,
const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
Array<float>& Dist, Array<float>& MultiScaleSmooth, Array<float>& NonLocalMean,
fillHalo<float>& fillFloat, const Domain& Dm, int nprocx, int level );
fillHalo<float>& fillFloat, const Domain& Dm, int nprocx, int level,
float threshold, float lamda, float sigsq, int depth);
// Remove regions that are likely noise by shrinking the volumes by dx,

22
common/Array.hpp
View File

@ -1064,15 +1064,15 @@ Array<TYPE, FUN> Array<TYPE, FUN>::reverseDim() const
template<class TYPE, class FUN>
Array<TYPE, FUN> Array<TYPE, FUN>::coarsen( const Array<TYPE, FUN> &filter ) const
{
auto S2 = size();
for ( size_t i = 0; i < S2.size(); i++ ) {
S2[i] /= filter.size( i );
if ( S2[i] * filter.size( i ) != size( i ) )
throw std::invalid_argument( "Array must be multiple of filter size" );
}
Array<TYPE, FUN> y( S2 );
if ( d_size.ndim() <= 3 )
throw std::logic_error( "Function programmed for more than 3 dimensions" );
auto S2 = size();
for ( size_t i = 0; i < S2.size(); i++ ) {
S2.resize( i, S2[i] / filter.size(i) );
if ( S2[i] * filter.size( i ) != size( i ) )
throw std::invalid_argument( "Array must be multiple of filter size" );
}
Array<TYPE, FUN> y( S2 );
if ( d_size.ndim() <= 3 )
throw std::logic_error( "Function programmed for more than 3 dimensions" );
const auto& Nh = filter.d_size;
for ( size_t k1 = 0; k1 < y.d_size[2]; k1++ ) {
for ( size_t j1 = 0; j1 < y.d_size[1]; j1++ ) {
@ -1096,8 +1096,8 @@ template<class TYPE, class FUN>
Array<TYPE, FUN> Array<TYPE, FUN>::coarsen(
const std::vector<size_t> &ratio, std::function<TYPE( const Array<TYPE, FUN> & )> filter ) const
{
if ( ratio.size() != d_size.ndim() )
throw std::logic_error( "ratio size does not match ndim" );
//if ( ratio.size() != d_size.ndim() )
// throw std::logic_error( "ratio size does not match ndim" );
auto S2 = size();
for ( size_t i = 0; i < S2.size(); i++ ) {
S2.resize( i, S2[i] / ratio[i] );

1
example/CMakeLists.txt
View File

@ -12,6 +12,7 @@ INSTALL_EXAMPLE( Poiseuille )
INSTALL_EXAMPLE( Juanes )
INSTALL_EXAMPLE( MicroModel )
INSTALL_EXAMPLE( Tiff )
INSTALL_EXAMPLE( uCT )
# Create unit tests for each example

32
example/uCT/eos-uCT.pbs Normal file
View File

@ -0,0 +1,32 @@
#!/bin/bash
#PBS -A GEO106
#PBS -N uCT
#PBS -j oe
##PBS -l walltime=02:30:00,nodes=20
#PBS -l walltime=00:10:00,nodes=20
##PBS -l gres=widow2%widow3
##PBS -q killable
#PBS -q debug
#cd /tmp/work/$USER
date
cd $PBS_O_WORKDIR
#LBPM_WIA_INSTALL_DIR=/lustre/atlas/proj-shared/geo106/build-eos-LBPM-WIA
LBPM_WIA_INSTALL_DIR=/ccs/proj/geo106/eos/LBPM-WIA
#echo "PBS_O_WORKDIR: `echo $PBS_O_WORKDIR`"
source $MODULESHOME/init/bash
module swap PrgEnv-intel PrgEnv-gnu
export LD_LIBRARY_PATH=${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}
# generate distance map from the binary image
# input files are ID.xxxxx
# output files are SignDist.xxxxx
NUMPROCS=320
aprun -n $NUMPROCS $LBPM_WIA_INSTALL_DIR/bin/lbpm_uCT_pp input.db
exit;

24
example/uCT/input.db Normal file
View File

@ -0,0 +1,24 @@
Domain {
nproc = 5, 8, 8 // Number of processors (Npx,Npy,Npz)
n = 240, 240, 240 // Size of local domain (Nx,Ny,Nz)
n_spheres = 1 // Number of spheres
L = 1, 1, 1 // Length of domain (x,y,z)
BC = 0 // Boundary condition type
}
uCT {
InputFile = "MG_dryrecon.volume" // this is the name of the data to segment
target = -160 // this is the target value of the phase to segment
background = 160 // this is the background value for all other phases
rough_cutoff = 0.0 // step 1 -- rough thresholding
lamda = 0.5 // step 2 -- apply some smoothing
nlm_sigsq = 0.25 // step 3 -- apply non-local means
nlm_depth = 3
center_x = 600
center_y = 960
center_z = 960
cylinder_radius = 700
}
Analysis {
}

2
sample_scripts/config_build_eos
View File

@ -36,7 +36,7 @@ cmake \
-D USE_CUDA=0 \
-D CUDA_FLAGS="-arch sm_35" \
-DBUILD_SHARED_LIBS=OFF \
-D USE_NETCDF=0 \
-D USE_NETCDF=1 \
-D NETCDF_DIRECTORY=$NETCDF_DIR \
-D HDF5_DIRECTORY=$HDF5_DIR \
-D CMAKE_SKIP_RPATH=true \

6
tests/CMakeLists.txt
View File

@ -67,13 +67,13 @@ ADD_LBPM_TEST( TestWriter )
IF ( USE_NETCDF )
ADD_LBPM_TEST_PARALLEL( TestNetcdf 8 )
ADD_LBPM_EXECUTABLE( lbpm_uCT_pp )
ADD_LBPM_EXECUTABLE( lbpm_uCT_maskfilter )
# ADD_LBPM_EXECUTABLE( lbpm_uCT_maskfilter )
ENDIF()
# Sample test that will run with 1, 2, and 4 processors, failing with 4 or more procs
ADD_LBPM_TEST_1_2_4( hello_world )
ADD_LBPM_TEST_1_2_4( TestColorBubble )
ADD_LBPM_TEST_1_2_4( TestColorSquareTube )
ADD_LBPM_TEST_1_2_4( TestColorBubble ../example/Bubble/input.db)
ADD_LBPM_TEST_1_2_4( TestColorSquareTube ../example/Bubble/input.db)
SET_TESTS_PROPERTIES( hello_world PROPERTIES ENVIRONMENT "MPICH_RDMA_ENABLED_CUDA=0")
IF ( USE_MPI )

76
tests/lbpm_segmented_pp.cpp
View File

@ -131,63 +131,35 @@ int main(int argc, char **argv)
// char fluidValue,solidValue;
int MAXTIME=1000;
int READ_FROM_BLOCK=0;
if (argc > 1){
MAXTIME=atoi(argv[1]);
}
if (argc > 2){
READ_FROM_BLOCK=atoi(argv[2]);
}
std::vector<char> solidValues;
std::vector<char> nwpValues;
std::string line;
if (rank==0){
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> nx;
domain >> ny;
domain >> nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
}
MPI_Barrier(comm);
// Computational domain
MPI_Bcast(&nx,1,MPI_INT,0,comm);
MPI_Bcast(&ny,1,MPI_INT,0,comm);
MPI_Bcast(&nz,1,MPI_INT,0,comm);
MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
//.................................................
MPI_Barrier(comm);
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
printf("Number of MPI ranks used: %i \n", nprocs);
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
fflush(stdout);
}
if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors");
}
char LocalRankString[8];
char LocalRankFilename[40];
string filename;
if (argc > 1) filename=argv[1];
else ERROR("No input database provided\n");
// read the input database
auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" );
// Read domain parameters
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
auto ReadValues = domain_db->getVector<char>( "ReadValues" );
auto WriteValues = domain_db->getVector<char>( "WriteValues" );
nx = size[0];
ny = size[1];
nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
int N = (nx+2)*(ny+2)*(nz+2);
std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
// std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
for (n=0; n<N; n++) Dm->id[n]=1;
Dm->CommInit();
std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );

228
tests/lbpm_serial_decomp.cpp
View File

@ -15,10 +15,10 @@
int main(int argc, char **argv)
{
int rank=0;
bool MULTINPUT=false;
/* bool MULTINPUT=false;
int NWP,SOLID,rank_offset;
SOLID=atoi(argv[1]);
@ -35,24 +35,62 @@ int main(int argc, char **argv)
MULTINPUT=true;
rank_offset=0;
}
*/
string filename;
if (argc > 1)
filename=argv[1];
else{
ERROR("lbpm_serial_decomp: no in put database provided \n");
}
int rank_offset=0;
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int64_t Nx,Ny,Nz;
int64_t i,j,k,n;
int BC=0;
char Filename[40];
int64_t xStart,yStart,zStart;
// char fluidValue,solidValue;
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocs, nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int64_t Nx,Ny,Nz;
int64_t i,j,k,n;
int BC=0;
int64_t xStart,yStart,zStart;
// char fluidValue,solidValue;
std::vector<char> solidValues;
std::vector<char> nwpValues;
std::string line;
xStart=yStart=zStart=0;
// read the input database
auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" );
if (rank==0){
// 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" );
auto ReadValues = domain_db->getVector<char>( "ReadValues" );
auto WriteValues = domain_db->getVector<char>( "WriteValues" );
nx = size[0];
ny = size[1];
nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
Nx = SIZE[0];
Ny = SIZE[1];
Nz = SIZE[2];
//Nx = nprocx*nx;
//Ny = nprocx*ny;
//Nz = nprocx*nz;
printf("Input media: %s\n",Filename.c_str());
printf("Relabeling %lu values\n",ReadValues.size());
for (int idx=0; idx<ReadValues.size(); idx++){
char oldvalue=ReadValues[idx];
char newvalue=WriteValues[idx];
printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue);
}
/* if (rank==0){
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
@ -61,7 +99,9 @@ int main(int argc, char **argv)
domain >> ny;
domain >> nz;
domain >> nspheres;
domain >> Lx;
domain >> Lx; printf("Domain decomposition completed successfully \n");
return 0;
domain >> Ly;
domain >> Lz;
@ -75,86 +115,100 @@ int main(int argc, char **argv)
image >> zStart;
}
nprocs=nprocx*nprocy*nprocz;
*/
nprocs=nprocx*nprocy*nprocz;
char *SegData = NULL;
// Rank=0 reads the entire segmented data and distributes to worker processes
if (rank==0){
printf("Dimensions of segmented image: %ld x %ld x %ld \n",Nx,Ny,Nz);
int64_t SIZE = Nx*Ny*Nz;
SegData = new char[SIZE];
FILE *SEGDAT = fopen(Filename,"rb");
if (SEGDAT==NULL) ERROR("Error reading segmented data");
size_t ReadSeg;
ReadSeg=fread(SegData,1,SIZE,SEGDAT);
if (ReadSeg != size_t(SIZE)) printf("lbpm_segmented_decomp: Error reading segmented data (rank=%i)\n",rank);
fclose(SEGDAT);
printf("Read segmented data from %s \n",Filename);
}
char *SegData = NULL;
// Rank=0 reads the entire segmented data and distributes to worker processes
if (rank==0){
printf("Dimensions of segmented image: %ld x %ld x %ld \n",Nx,Ny,Nz);
int64_t SIZE = Nx*Ny*Nz;
SegData = new char[SIZE];
FILE *SEGDAT = fopen(Filename.c_str(),"rb");
if (SEGDAT==NULL) ERROR("Error reading segmented data");
size_t ReadSeg;
ReadSeg=fread(SegData,1,SIZE,SEGDAT);
if (ReadSeg != size_t(SIZE)) printf("lbpm_segmented_decomp: Error reading segmented data (rank=%i)\n",rank);
fclose(SEGDAT);
printf("Read segmented data from %s \n",Filename.c_str());
}
// Get the rank info
int64_t N = (nx+2)*(ny+2)*(nz+2);
// 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 - (Nz - zStart))/2;
if (z_transition_size < 0) z_transition_size=0;
// number of sites to use for periodic boundary condition transition zone
int64_t z_transition_size = (nprocz*nz - (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)];
char LocalRankFilename[40];
char *loc_id;
loc_id = new char [(nx+2)*(ny+2)*(nz+2)];
// 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);
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 (z<zStart) z=zStart;
if (!(z<Nz)) z=Nz-1;
int64_t nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
int64_t nglobal = z*Nx*Ny+y*Nx+x;
loc_id[nlocal] = SegData[nglobal];
}
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 (z<zStart) z=zStart;
if (!(z<Nz)) z=Nz-1;
int64_t nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
int64_t nglobal = z*Nx*Ny+y*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;;
if (loc_id[n]==char(SOLID)) loc_id[n] = 0;
else if (loc_id[n]==char(NWP)) loc_id[n] = 1;
else loc_id[n] = 2;
}
}
}
// 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);
}
// 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++){
char oldvalue=ReadValues[idx];
char newvalue=WriteValues[idx];
if (locval == oldvalue){
loc_id[n] = newvalue;
LabelCount[idx]++;
idx = ReadValues.size();
}
}
//if (loc_id[n]==char(SOLID)) loc_id[n] = 0;
//else if (loc_id[n]==char(NWP)) loc_id[n] = 1;
//else loc_id[n] = 2;
}
}
}
// 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);
}
}
}
printf("Domain decomposition completed successfully \n");
return 0;
}
for (int idx=0; idx<ReadValues.size(); idx++){
char label=ReadValues[idx];
int count=LabelCount[idx];
printf("Label=%d, Count=%d \n",label,count);
}
}

630
tests/lbpm_uCT_pp.cpp
View File

@ -20,13 +20,12 @@
#include "IO/Writer.h"
#include "IO/netcdf.h"
#include "analysis/analysis.h"
#include "analysis/eikonal.h"
#include "analysis/filters.h"
#include "analysis/uCT.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
int main(int argc, char **argv)
{
@ -36,238 +35,303 @@ int main(int argc, char **argv)
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
Utilities::setErrorHandlers();
PROFILE_START("Main");
{
Utilities::setErrorHandlers();
PROFILE_START("Main");
//std::vector<std::string> filenames;
if ( argc<2 ) {
if ( rank == 0 )
printf("At least one filename must be specified\n");
return 1;
}
std::string filename = std::string(argv[1]);
if ( rank == 0 )
printf("Input data file: %s\n",filename.c_str());
//std::vector<std::string> filenames;
if ( argc<2 ) {
if ( rank == 0 ){
printf("At least one filename must be specified\n");
}
return 1;
}
std::string filename = std::string(argv[1]);
if ( rank == 0 ){
printf("Input data file: %s\n",filename.c_str());
}
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
read_domain( rank, nprocs, comm, nprocx, nprocy, nprocz, nx, ny, nz, nspheres, Lx, Ly, Lz );
int BC=0;
auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" );
auto uct_db = db->getDatabase( "uCT" );
auto analysis_db = db->getDatabase( "Analysis" );
// Read domain values
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
int BoundaryCondition = domain_db->getScalar<int>( "BC" );
int nx = size[0];
int ny = size[1];
int nz = size[2];
double Lx = L[0];
double Ly = L[1];
double Lz = L[2];
int nprocx = nproc[0];
int nprocy = nproc[1];
int nprocz = nproc[2];
auto InputFile=uct_db->getScalar<std::string>( "InputFile" );
auto target=uct_db->getScalar<int>("target");
auto background=uct_db->getScalar<int>("background");
auto rough_cutoff=uct_db->getScalar<float>( "rough_cutoff" );
auto lamda=uct_db->getScalar<float>( "lamda" );
auto nlm_sigsq=uct_db->getScalar<float>( "nlm_sigsq" );
auto nlm_depth=uct_db->getScalar<int>( "nlm_depth" );
auto cx=uct_db->getScalar<int>( "center_x" );
auto cy=uct_db->getScalar<int>( "center_y" );
auto cz=uct_db->getScalar<int>( "center_z" );
auto CylRad=uct_db->getScalar<float>( "cylinder_radius" );
//.......................................................................
// Reading the domain information file
//.......................................................................
// std::shared_ptr<Domain> Dm ();
//for (int i=0; i<Dm->Nx*Dm->Ny*Dm->Nz; i++) Dm->id[i] = 1;
//Dm->CommInit();
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
printf("Number of MPI ranks used: %i \n", nprocs);
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
}
if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors");
}
// Determine the maximum number of levels for the desired coarsen ratio
int ratio[3] = {2,2,2};
//std::vector<size_t> ratio = {4,4,4};
// need to set up databases for each level of the mesh
std:vector<Database> multidomain_db;
std::vector<int> Nx(1,nx), Ny(1,ny), Nz(1,nz);
while ( Nx.back()%ratio[0]==0 && Nx.back()>8 &&
Ny.back()%ratio[1]==0 && Ny.back()>8 &&
Nz.back()%ratio[2]==0 && Nz.back()>8 )
{
Nx.push_back( Nx.back()/ratio[0] );
Ny.push_back( Ny.back()/ratio[1] );
Nz.push_back( Nz.back()/ratio[2] );
// clone the domain and create coarse version based on Nx,Ny,Nz
//multidomain_db.push_back();
}
int N_levels = Nx.size();
// Initialize the domain
std::vector<std::shared_ptr<Domain>> Dm(N_levels);
for (int i=0; i<N_levels; i++) {
// This line is no good -- will create identical Domain structures instead of
// Need a way to define a coarse structure for the coarse domain (see above)
Dm[i].reset( new Domain(domain_db, comm) );
int N = (Nx[i]+2)*(Ny[i]+2)*(Nz[i]+2);
for (int n=0; n<N; n++){
Dm[i]->id[n] = 1;
}
Dm[i]->CommInit();
}
// array containing a distance mask
Array<float> MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2);
// Create the level data
std::vector<Array<char>> ID(N_levels);
std::vector<Array<float>> LOCVOL(N_levels);
std::vector<Array<float>> Dist(N_levels);
std::vector<Array<float>> MultiScaleSmooth(N_levels);
std::vector<Array<float>> Mean(N_levels);
std::vector<Array<float>> NonLocalMean(N_levels);
std::vector<std::shared_ptr<fillHalo<double>>> fillDouble(N_levels);
std::vector<std::shared_ptr<fillHalo<float>>> fillFloat(N_levels);
std::vector<std::shared_ptr<fillHalo<char>>> fillChar(N_levels);
for (int i=0; i<N_levels; i++) {
ID[i] = Array<char>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
LOCVOL[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
Dist[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
MultiScaleSmooth[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
Mean[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
NonLocalMean[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
ID[i].fill(0);
LOCVOL[i].fill(0);
Dist[i].fill(0);
MultiScaleSmooth[i].fill(0);
Mean[i].fill(0);
NonLocalMean[i].fill(0);
fillDouble[i].reset(new fillHalo<double>(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) );
fillFloat[i].reset(new fillHalo<float>(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) );
fillChar[i].reset(new fillHalo<char>(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) );
}
// Read the subvolume of interest on each processor
PROFILE_START("ReadVolume");
int fid = netcdf::open(InputFile,netcdf::READ);
std::string varname("VOLUME");
netcdf::VariableType type = netcdf::getVarType( fid, varname );
std::vector<size_t> dim = netcdf::getVarDim( fid, varname );
if ( rank == 0 ) {
printf("Reading %s (%s)\n",varname.c_str(),netcdf::VariableTypeName(type).c_str());
printf(" dims = %i x %i x %i \n",int(dim[0]),int(dim[1]),int(dim[2]));
}
{
RankInfoStruct info( rank, nprocx, nprocy, nprocz );
int x = info.ix*nx;
int y = info.jy*ny;
int z = info.kz*nz;
// Read the local data
Array<short> VOLUME = netcdf::getVar<short>( fid, varname, {x,y,z}, {nx,ny,nz}, {1,1,1} );
// Copy the data and fill the halos
LOCVOL[0].fill(0);
fillFloat[0]->copy( VOLUME, LOCVOL[0] );
fillFloat[0]->fill( LOCVOL[0] );
}
netcdf::close( fid );
MPI_Barrier(comm);
PROFILE_STOP("ReadVolume");
if (rank==0) printf("Read complete\n");
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
printf("Number of MPI ranks used: %i \n", nprocs);
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
}
if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors");
}
// Filter the original data
filter_src( *Dm[0], LOCVOL[0] );
// Determine the maximum number of levels for the desired coarsen ratio
int ratio[3] = {4,4,4};
std::vector<int> Nx(1,nx), Ny(1,ny), Nz(1,nz);
while ( Nx.back()%ratio[0]==0 && Nx.back()>8 &&
Ny.back()%ratio[1]==0 && Ny.back()>8 &&
Nz.back()%ratio[2]==0 && Nz.back()>8 )
{
Nx.push_back( Nx.back()/ratio[0] );
Ny.push_back( Ny.back()/ratio[1] );
Nz.push_back( Nz.back()/ratio[2] );
}
int N_levels = Nx.size();
// Set up the mask to be distance to cylinder (crop outside cylinder)
// float CylRad=900;
for (int k=0;k<Nz[0]+2;k++) {
for (int j=0;j<Ny[0]+2;j++) {
for (int i=0;i<Nx[0]+2;i++) {
int iproc = Dm[0]->iproc();
int jproc = Dm[0]->jproc();
int kproc = Dm[0]->kproc();
// Initialize the domain
std::vector<std::shared_ptr<Domain>> Dm(N_levels);
for (int i=0; i<N_levels; i++) {
Dm[i].reset( new Domain(Nx[i],Ny[i],Nz[i],rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC) );
int N = (Nx[i]+2)*(Ny[i]+2)*(Nz[i]+2);
for (int n=0; n<N; n++)
Dm[i]->id[n] = 1;
Dm[i]->CommInit(comm);
}
int x=iproc*Nx[0]+i-1;
int y=jproc*Ny[0]+j-1;
int z=kproc*Nz[0]+k-1;
// array containing a distance mask
Array<float> MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2);
// Create the level data
std::vector<Array<char>> ID(N_levels);
std::vector<Array<float>> LOCVOL(N_levels);
std::vector<Array<float>> Dist(N_levels);
std::vector<Array<float>> MultiScaleSmooth(N_levels);
std::vector<Array<float>> Mean(N_levels);
std::vector<Array<float>> NonLocalMean(N_levels);
std::vector<std::shared_ptr<fillHalo<double>>> fillDouble(N_levels);
std::vector<std::shared_ptr<fillHalo<float>>> fillFloat(N_levels);
std::vector<std::shared_ptr<fillHalo<char>>> fillChar(N_levels);
for (int i=0; i<N_levels; i++) {
ID[i] = Array<char>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
LOCVOL[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
Dist[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
MultiScaleSmooth[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
Mean[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
NonLocalMean[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
ID[i].fill(0);
LOCVOL[i].fill(0);
Dist[i].fill(0);
MultiScaleSmooth[i].fill(0);
Mean[i].fill(0);
NonLocalMean[i].fill(0);
fillDouble[i].reset(new fillHalo<double>(Dm[i]->Comm,Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],1,1,1,0,1) );
fillFloat[i].reset(new fillHalo<float>(Dm[i]->Comm,Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],1,1,1,0,1) );
fillChar[i].reset(new fillHalo<char>(Dm[i]->Comm,Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],1,1,1,0,1) );
}
//int cx = 0.5*nprocx*Nx[0];
//int cy = 0.5*nprocy*Ny[0];
//int cz = 0.5*nprocz*Nz[0];
// distance from the center line
MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) );
// Read the subvolume of interest on each processor
PROFILE_START("ReadVolume");
int fid = netcdf::open(filename,netcdf::READ);
std::string varname("VOLUME");
netcdf::VariableType type = netcdf::getVarType( fid, varname );
std::vector<size_t> dim = netcdf::getVarDim( fid, varname );
if ( rank == 0 ) {
printf("Reading %s (%s)\n",varname.c_str(),netcdf::VariableTypeName(type).c_str());
printf(" dims = %i x %i x %i \n",int(dim[0]),int(dim[1]),int(dim[2]));
}
{
RankInfoStruct info( rank, nprocx, nprocy, nprocz );
int x = info.ix*nx;
int y = info.jy*ny;
int z = info.kz*nz;
// Read the local data
Array<short> VOLUME = netcdf::getVar<short>( fid, varname, {x,y,z}, {nx,ny,nz}, {1,1,1} );
// Copy the data and fill the halos
LOCVOL[0].fill(0);
fillFloat[0]->copy( VOLUME, LOCVOL[0] );
fillFloat[0]->fill( LOCVOL[0] );
}
netcdf::close( fid );
MPI_Barrier(comm);
PROFILE_STOP("ReadVolume");
if (rank==0) printf("Read complete\n");
// Filter the original data
filter_src( *Dm[0], LOCVOL[0] );
// Set up the mask to be distance to cylinder (crop outside cylinder)
float CylRad=900;
for (int k=0;k<Nz[0]+2;k++) {
for (int j=0;j<Ny[0]+2;j++) {
for (int i=0;i<Nx[0]+2;i++) {
int iproc = Dm[0]->iproc;
int jproc = Dm[0]->jproc;
int kproc = Dm[0]->kproc;
int x=iproc*Nx[0]+i-1;
int y=jproc*Ny[0]+j-1;
int z=kproc*Nz[0]+k-1;
int cx = 0.5*nprocx*Nx[0];
int cy = 0.5*nprocy*Ny[0];
int cz = 0.5*nprocz*Nz[0];
// distance from the center line
MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) );
}
}
}
}
// Compute the means for the high/low regions
// (should use automated mixture model to approximate histograms)
float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) );
float mean_plus=0;
float mean_minus=0;
int count_plus=0;
int count_minus=0;
for (int k=1;k<Nz[0]+1;k++) {
for (int j=1;j<Ny[0]+1;j++) {
for (int i=1;i<Nx[0]+1;i++) {
if (MASK(i,j,k) > 0.f ){
auto tmp = LOCVOL[0](i,j,k);
if ( tmp > THRESHOLD ) {
mean_plus += tmp;
count_plus++;
} else if ( tmp < -THRESHOLD ) {
mean_minus += tmp;
count_minus++;
// Compute the means for the high/low regions
// (should use automated mixture model to approximate histograms)
//float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) );
float THRESHOLD=0.5*float(target+background);
float mean_plus=0;
float mean_minus=0;
int count_plus=0;
int count_minus=0;
for (int k=1;k<Nz[0]+1;k++) {
for (int j=1;j<Ny[0]+1;j++) {
for (int i=1;i<Nx[0]+1;i++) {
if (MASK(i,j,k) > 0.f ){
auto tmp = LOCVOL[0](i,j,k);
/* if ((tmp-background)*(tmp-target) > 0){
// direction to background / target is the same
if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background
else tmp=target; // tmp closer to target
}
*/
if ( tmp > THRESHOLD ) {
mean_plus += tmp;
count_plus++;
} else if ( tmp < -THRESHOLD ) {
mean_minus += tmp;
count_minus++;
}
}
}
}
}
mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus );
mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus );
if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus);
MPI_Barrier(comm);
// Scale the source data to +-1.0
for (size_t i=0; i<LOCVOL[0].length(); i++) {
if (MASK(i) < 0.f){
LOCVOL[0](i) = 1.0;
}
else if ( LOCVOL[0](i) >= 0 ) {
LOCVOL[0](i) /= mean_plus;
LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f );
} else {
LOCVOL[0](i) /= -mean_minus;
LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f );
}
}
// Fill the source data for the coarse meshes
PROFILE_START("CoarsenMesh");
for (int i=1; i<N_levels; i++) {
Array<float> filter(ratio[0],ratio[1],ratio[2]);
filter.fill(1.0f/filter.length());
Array<float> tmp(Nx[i-1],Ny[i-1],Nz[i-1]);
fillFloat[i-1]->copy( LOCVOL[i-1], tmp );
Array<float> coarse = tmp.coarsen( filter );
fillFloat[i]->copy( coarse, LOCVOL[i] );
fillFloat[i]->fill( LOCVOL[i] );
if (rank==0){
printf("Coarsen level %i \n",i);
printf(" Nx=%i, Ny=%i, Nz=%i \n",int(tmp.size(0)),int(tmp.size(1)),int(tmp.size(2)) );
printf(" filter_x=%i, filter_y=%i, filter_z=%i \n",int(filter.size(0)),int(filter.size(1)),int(filter.size(2)) );
printf(" ratio= %i,%i,%i \n",int(ratio[0]),int(ratio[1]),int(ratio[2]) );
}
MPI_Barrier(comm);
}
PROFILE_STOP("CoarsenMesh");
// Initialize the coarse level
PROFILE_START("Solve full mesh");
if (rank==0)
printf("Initialize full mesh\n");
solve( LOCVOL[0], Mean[0], ID[0], Dist[0], MultiScaleSmooth[0],
NonLocalMean[0], *fillFloat[0], *Dm[0], nprocx,
rough_cutoff, lamda, nlm_sigsq, nlm_depth);
PROFILE_STOP("Solve full mesh");
MPI_Barrier(comm);
/* // Initialize the coarse level
PROFILE_START("Solve coarse mesh");
if (rank==0)
printf("Initialize coarse mesh\n");
solve( LOCVOL.back(), Mean.back(), ID.back(), Dist.back(), MultiScaleSmooth.back(),
NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx );
PROFILE_STOP("Solve coarse mesh");
MPI_Barrier(comm);
// Refine the solution
PROFILE_START("Refine distance");
if (rank==0)
printf("Refine mesh\n");
for (int i=int(Nx.size())-2; i>=0; i--) {
if (rank==0)
printf(" Refining to level %i\n",int(i));
refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i],
NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i );
}
mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus );
mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus );
if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus);
PROFILE_STOP("Refine distance");
MPI_Barrier(comm);
// Perform a final filter
PROFILE_START("Filtering final domains");
if (rank==0)
printf("Filtering final domains\n");
Array<float> filter_Mean, filter_Dist1, filter_Dist2;
filter_final( ID[0], Dist[0], *fillFloat[0], *Dm[0], filter_Mean, filter_Dist1, filter_Dist2 );
PROFILE_STOP("Filtering final domains");
*/
//removeDisconnected( ID[0], *Dm[0] );
// Scale the source data to +-1.0
for (size_t i=0; i<LOCVOL[0].length(); i++) {
if (MASK(i) < 0.f){
LOCVOL[0](i) = 1.0;
}
else if ( LOCVOL[0](i) >= 0 ) {
LOCVOL[0](i) /= mean_plus;
LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f );
} else {
LOCVOL[0](i) /= -mean_minus;
LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f );
}
}
// Fill the source data for the coarse meshes
PROFILE_START("CoarsenMesh");
for (int i=1; i<N_levels; i++) {
Array<float> filter(ratio[0],ratio[1],ratio[2]);
filter.fill(1.0f/filter.length());
Array<float> tmp(Nx[i-1],Ny[i-1],Nz[i-1]);
fillFloat[i-1]->copy( LOCVOL[i-1], tmp );
Array<float> coarse = tmp.coarsen( filter );
fillFloat[i]->copy( coarse, LOCVOL[i] );
fillFloat[i]->fill( LOCVOL[i] );
}
PROFILE_STOP("CoarsenMesh");
// Initialize the coarse level
PROFILE_START("Solve coarse mesh");
if (rank==0)
printf("Initialize coarse mesh\n");
solve( LOCVOL.back(), Mean.back(), ID.back(), Dist.back(), MultiScaleSmooth.back(),
NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx );
PROFILE_STOP("Solve coarse mesh");
// Refine the solution
PROFILE_START("Refine distance");
if (rank==0)
printf("Refine mesh\n");
for (int i=int(Nx.size())-2; i>=0; i--) {
if (rank==0)
printf(" Refining to level %i\n",int(i));
refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i],
NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i );
}
PROFILE_STOP("Refine distance");
// Perform a final filter
PROFILE_START("Filtering final domains");
if (rank==0)
printf("Filtering final domains\n");
Array<float> filter_Mean, filter_Dist1, filter_Dist2;
filter_final( ID[0], Dist[0], *fillFloat[0], *Dm[0], filter_Mean, filter_Dist1, filter_Dist2 );
PROFILE_STOP("Filtering final domains");
//removeDisconnected( ID[0], *Dm[0] );
// Write the distance function to a netcdf file
/*
// Write the distance function to a netcdf file
const char* netcdf_filename = "Distance.nc";
{
RankInfoStruct info( rank, nprocx, nprocy, nprocz );
@ -279,33 +343,59 @@ int main(int argc, char **argv)
netcdf::write( fid, "Distance", dims, data, info );
netcdf::close( fid );
}
*/
{
// Write the results to visit
if (rank==0) printf("Setting up visualization structure \n");
// std::vector<IO::MeshDataStruct> meshData(N_levels);
std::vector<IO::MeshDataStruct> meshData(1);
// for (size_t i=0; i<Nx.size(); i++) {
// Mesh
meshData[0].meshName = "image";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm[0]->rank_info,Nx[0],Ny[0],Nz[0],Lx,Ly,Lz) );
// Source data
std::shared_ptr<IO::Variable> OrigData( new IO::Variable() );
OrigData->name = "Source Data";
OrigData->type = IO::VariableType::VolumeVariable;
OrigData->dim = 1;
OrigData->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(OrigData);
fillDouble[0]->copy( LOCVOL[0], OrigData->data );
// Non-Local Mean
std::shared_ptr<IO::Variable> NonLocMean( new IO::Variable() );
NonLocMean->name = "Non-Local Mean";
NonLocMean->type = IO::VariableType::VolumeVariable;
NonLocMean->dim = 1;
NonLocMean->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(NonLocMean);
fillDouble[0]->copy( NonLocalMean[0], NonLocMean->data );
std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
SegData->name = "Segmented Data";
SegData->type = IO::VariableType::VolumeVariable;
SegData->dim = 1;
SegData->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(SegData);
fillDouble[0]->copy( ID[0], SegData->data );
// Signed Distance
std::shared_ptr<IO::Variable> DistData( new IO::Variable() );
DistData->name = "Signed Distance";
DistData->type = IO::VariableType::VolumeVariable;
DistData->dim = 1;
DistData->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(DistData);
fillDouble[0]->copy( Dist[0], DistData->data );
// Smoothed Data
std::shared_ptr<IO::Variable> SmoothData( new IO::Variable() );
SmoothData->name = "Smoothed Data";
SmoothData->type = IO::VariableType::VolumeVariable;
SmoothData->dim = 1;
SmoothData->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(SmoothData);
fillDouble[0]->copy( MultiScaleSmooth[0], SmoothData->data );
// Write the results to visit
if (rank==0) printf("Writing output \n");
std::vector<IO::MeshDataStruct> meshData(N_levels);
for (size_t i=0; i<Nx.size(); i++) {
// Mesh
meshData[i].meshName = "Level " + std::to_string(i+1);
meshData[i].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],Lx,Ly,Lz) );
// Source data
std::shared_ptr<IO::Variable> OrigData( new IO::Variable() );
OrigData->name = "Source Data";
OrigData->type = IO::VariableType::VolumeVariable;
OrigData->dim = 1;
OrigData->data.resize(Nx[i],Ny[i],Nz[i]);
meshData[i].vars.push_back(OrigData);
fillDouble[i]->copy( LOCVOL[i], OrigData->data );
// Non-Local Mean
std::shared_ptr<IO::Variable> NonLocMean( new IO::Variable() );
NonLocMean->name = "Non-Local Mean";
NonLocMean->type = IO::VariableType::VolumeVariable;
NonLocMean->dim = 1;
NonLocMean->data.resize(Nx[i],Ny[i],Nz[i]);
meshData[i].vars.push_back(NonLocMean);
fillDouble[i]->copy( NonLocalMean[i], NonLocMean->data );
// Segmented Data
std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
/*// Segmented Data
std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
SegData->name = "Segmented Data";
SegData->type = IO::VariableType::VolumeVariable;
SegData->dim = 1;
@ -328,8 +418,9 @@ int main(int argc, char **argv)
SmoothData->data.resize(Nx[i],Ny[i],Nz[i]);
meshData[i].vars.push_back(SmoothData);
fillDouble[i]->copy( MultiScaleSmooth[i], SmoothData->data );
}
#if 0
#if 0
std::shared_ptr<IO::Variable> filter_Mean_var( new IO::Variable() );
filter_Mean_var->name = "Mean";
filter_Mean_var->type = IO::VariableType::VolumeVariable;
@ -352,37 +443,40 @@ int main(int argc, char **argv)
meshData[0].vars.push_back(filter_Dist2_var);
fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data );
#endif
// Write visulization data
IO::writeData( 0, meshData, comm );
if (rank==0) printf("Finished. \n");
/* 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 (Dm.id[n]==char(SOLID)) Dm.id[n] = 0;
else if (Dm.id[n]==char(NWP)) Dm.id[n] = 1;
else Dm.id[n] = 2;
*/
MPI_Barrier(comm);
if (rank==0) printf("Writing output \n");
// Write visulization data
IO::writeData( 0, meshData, comm );
if (rank==0) printf("Finished. \n");
}
// Compute the Minkowski functionals
MPI_Barrier(comm);
std::shared_ptr<Minkowski> Averages(new Minkowski(Dm[0]));
if (rank==0) printf("Initializing the system: Nx=%i, Ny=%i, Nz=%i \n",Averages->Nx,Averages->Ny,Averages->Nz);
for ( int k=1;k<Averages->Nx-1;k++){
for ( int j=1;j<Averages->Ny-1;j++){
for ( int i=1;i<Averages->Nx-1;i++){
int n = k*(Averages->Nx)*(Averages->Ny)+j*(Averages->Nx)+i;
double distance=Dist[0](i,j,k);
Averages->SDn(i,j,k) = distance;
if (distance > 0)
Averages->Dm->id[n] = 0;
else
Averages->Dm->id[n] = 1;
}
}
}
if (rank==0) printf("Computing Minkowski functionals \n");
Averages->Initialize();
Averages->UpdateMeshValues();
Averages->ComputeLocal();
Averages->Reduce();
Averages->PrintAll();
}
if (rank==0) printf("Domain set \n");
// Write the local volume files
char LocalRankString[8];
char LocalRankFilename[40];
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"Seg.%s",LocalRankString);
FILE * SEG;
SEG=fopen(LocalRankFilename,"wb");
fwrite(LOCVOL.get(),4,N,SEG);
fclose(SEG);
*/
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_uCT_pp",true);
PROFILE_SAVE("lbpm_uCT_pp",true);
MPI_Barrier(comm);
MPI_Finalize();
return 0;