Updating lbpm_uCT_pp.cpp

IO/Writer.h

@@ -24,7 +24,7 @@ namespace IO {
  *                          silo - Silo
  * @param[in] append        Append any existing data (default is false)
  */
-void initialize( const std::string& path="", const std::string& format="new", bool append=false );
+void initialize( const std::string& path="", const std::string& format="silo", bool append=false );
 
 
 /*!

analysis/runAnalysis.cpp

@@ -23,8 +23,8 @@ AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs)
 }
 bool matches( AnalysisType x, AnalysisType y )
 {
-	return static_cast<std::underlying_type<AnalysisType>::type>(x) &
+    return ( static_cast<std::underlying_type<AnalysisType>::type>(x) &
-			static_cast<std::underlying_type<AnalysisType>::type>(y) != 0;
+             static_cast<std::underlying_type<AnalysisType>::type>(y) ) != 0;
 }
 
 
@@ -40,7 +40,7 @@ class WriteRestartWorkItem: public ThreadPool::WorkItemRet<void>
 {
 public:
     WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cDen_, std::shared_ptr<double> cfq_, int N_ ):
-		filename(filename_), cDen(cDen_), cfq(cfq_), N(N_) {}
+        filename(filename_), cfq(cfq_), cDen(cDen_), N(N_) {}
     virtual void run() {
         PROFILE_START("Save Checkpoint",1);
         double value;
@@ -303,10 +303,16 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
             d_regular ( Regular),
             d_rank_info( rank_info ),
             d_Map( Map ),
-    		d_ScaLBL_Comm( ScaLBL_Comm),
+            d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1),
-    		d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1)
+            d_ScaLBL_Comm( ScaLBL_Comm)
 {
 
+    // Ids of work items to use for dependencies
+    ThreadPool::thread_id_t d_wait_blobID;
+    ThreadPool::thread_id_t d_wait_analysis;
+    ThreadPool::thread_id_t d_wait_vis;
+    ThreadPool::thread_id_t d_wait_restart;
+
     char rankString[20];
     sprintf(rankString,"%05d",Dm->rank());
     d_N[0] = Dm->Nx;

common/Array.hpp

@@ -1002,7 +1002,7 @@ Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::coarsen(
             throw std::invalid_argument( "Array must be multiple of filter size" );
     }
     Array<TYPE, FUN, Allocator> y( S2 );
-    if ( d_size.ndim() <= 3 )
+    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++ ) {

tests/lbpm_uCT_pp.cpp

@@ -79,17 +79,12 @@ int main(int argc, char **argv)
         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 center       = uct_db->getVector<int>( "center" );
-		auto cy=uct_db->getScalar<int>( "center_y" );
-		auto cz=uct_db->getScalar<int>( "center_z" );
         auto CylRad       = uct_db->getScalar<float>( "cylinder_radius" );
-
+        auto maxLevels    = uct_db->getScalar<int>( "max_levels" );
-		//.......................................................................
+        std::vector<int> offset( 3, 0 );
-		// Reading the domain information file
+        if ( uct_db->keyExists( "offset" ) )
-		//.......................................................................
+            offset = uct_db->getVector<int>( "offset" );
-		//    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 ) {
@@ -105,18 +100,20 @@ int main(int argc, char **argv)
         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<std::shared_ptr<Database>> multidomain_db(1,domain_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 )
+                Nz.back()%ratio[2]==0 && Nz.back()>8 &&
+                (int) Nx.size() < maxLevels )
         {
             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();
+            auto db2 = domain_db->cloneDatabase();
+            db2->putVector<int>( "n", { Nx.back(), Ny.back(), Nz.back() } );
+            multidomain_db.push_back(db2);
         }
         int N_levels = Nx.size();
 
@@ -125,7 +122,7 @@ int main(int argc, char **argv)
         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) );
+            Dm[i].reset( new Domain(multidomain_db[i], comm) );
             int N = (Nx[i]+2)*(Ny[i]+2)*(Nz[i]+2);
             for (int n=0; n<N; n++){
                 Dm[i]->id[n] = 1;
@@ -135,6 +132,7 @@ int main(int argc, char **argv)
 
         // array containing a distance mask
         Array<float> MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2);
+        MASK.fill(0);
 
         // Create the level data
         std::vector<Array<char>>  ID(N_levels);
@@ -168,18 +166,17 @@ int main(int argc, char **argv)
         PROFILE_START("ReadVolume");
         int fid = netcdf::open(InputFile,netcdf::READ);
         std::string varname("VOLUME");
-		netcdf::VariableType type = netcdf::getVarType( fid, varname );
+        auto type = netcdf::getVarType( fid, varname );
-		std::vector<size_t> dim = netcdf::getVarDim( fid, varname );
+        auto 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
+            int x = Dm[0]->iproc()*nx + offset[0];
+            int y = Dm[0]->jproc()*ny + offset[1];
+            int z = Dm[0]->kproc()*nz + offset[2];
             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);
@@ -196,25 +193,17 @@ int main(int argc, char **argv)
         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 x=Dm[0]->iproc()*Nx[0]+i-1;
-					int jproc = Dm[0]->jproc();
+                    int y=Dm[0]->jproc()*Ny[0]+j-1;
-					int kproc = Dm[0]->kproc();
+                    int z=Dm[0]->kproc()*Nz[0]+k-1;
-
+                    int cx = center[0] - offset[0];
-					int x=iproc*Nx[0]+i-1;
+                    int cy = center[1] - offset[1];
-					int y=jproc*Ny[0]+j-1;
+                    int cz = center[2] - offset[2];
-					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)) );
-
                 }
             }
         }
@@ -222,9 +211,9 @@ int main(int argc, char **argv)
         // 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);
+        double THRESHOLD=0.5*(target+background);
-		float mean_plus=0;
+        double mean_plus=0;
-		float mean_minus=0;
+        double mean_minus=0;
         int count_plus=0;
         int count_minus=0;
         for (int k=1;k<Nz[0]+1;k++) {
@@ -249,6 +238,7 @@ int main(int argc, char **argv)
                 }
             }
         }
+        ASSERT( count_plus > 0 && count_minus > 0 );
         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);
@@ -269,6 +259,7 @@ int main(int argc, char **argv)
             }
         }
 
+
         // Fill the source data for the coarse meshes
         PROFILE_START("CoarsenMesh");
         for (int i=1; i<N_levels; i++) {
@@ -290,21 +281,12 @@ int main(int argc, char **argv)
         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");
+            printf("Initialize full mesh\n");
         solve( LOCVOL.back(), Mean.back(), ID.back(), Dist.back(), MultiScaleSmooth.back(),
-			NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx );
+                NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx, 
+                rough_cutoff, lamda, nlm_sigsq, nlm_depth);
         PROFILE_STOP("Solve coarse mesh");
         MPI_Barrier(comm);
 
@@ -312,11 +294,12 @@ int main(int argc, char **argv)
         PROFILE_START("Refine distance");
         if (rank==0)
             printf("Refine mesh\n");
-	for (int i=int(Nx.size())-2; i>=0; i--) {
+        for (int i=N_levels-2; i>=0; i--) {
             if (rank==0)
-			printf("   Refining to level %i\n",int(i));
+                printf("   Refining to level %i\n",i);
             refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i],
-				NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i );
+                    NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i, 
+                rough_cutoff, lamda, nlm_sigsq, nlm_depth);
         }
         PROFILE_STOP("Refine distance");
         MPI_Barrier(comm);    
@@ -328,12 +311,10 @@ int main(int argc, char **argv)
         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
-    const char* netcdf_filename = "Distance.nc";
+        /* const char* netcdf_filename = "Distance.nc";
         {
             RankInfoStruct info( rank, nprocx, nprocy, nprocz );
             std::vector<int> dim = { Nx[0]*nprocx, Ny[0]*nprocy, Nz[0]*nprocz };
@@ -343,77 +324,50 @@ int main(int argc, char **argv)
             fillFloat[0]->copy( Dist[0], data );
             netcdf::write( fid, "Distance", dims, data, info );
             netcdf::close( fid );
-    }
+        } */
-		 */
 
-		{
         // Write the results
         if (rank==0) printf("Setting up visualization structure \n");
-			//	std::vector<IO::MeshDataStruct> meshData(N_levels);
+        std::vector<IO::MeshDataStruct> meshData(N_levels);
-			std::vector<IO::MeshDataStruct> meshData(1);
+        for (size_t i=0; i<Nx.size(); i++) {
-			//    for (size_t i=0; i<Nx.size(); i++) {
             // Mesh
-			meshData[0].meshName = "image";
+            meshData[i].meshName = "image_" + std::to_string( i );
-			meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm[0]->rank_info,Nx[0],Ny[0],Nz[0],Lx,Ly,Lz) );
+            meshData[i].mesh = std::make_shared<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() );
+            auto OrigData = std::make_shared<IO::Variable>();
-			OrigData->name = "Source Data";
+            OrigData->name = "Source_Data_" + std::to_string( i );
             OrigData->type = IO::VariableType::VolumeVariable;
             OrigData->dim = 1;
-			OrigData->data.resize(Nx[0],Ny[0],Nz[0]);
+            OrigData->data.resize(Nx[i],Ny[i],Nz[i]);
-			meshData[0].vars.push_back(OrigData);
+            meshData[i].vars.push_back(OrigData);
-			fillDouble[0]->copy( LOCVOL[0], OrigData->data );
+            fillDouble[i]->copy( LOCVOL[i], OrigData->data );
             // Non-Local Mean
-			std::shared_ptr<IO::Variable> NonLocMean( new IO::Variable() );
+            auto NonLocMean = std::make_shared<IO::Variable>();
-			NonLocMean->name = "Non-Local Mean";
+            NonLocMean->name = "NonLocal_Mean_" + std::to_string( i );
             NonLocMean->type = IO::VariableType::VolumeVariable;
             NonLocMean->dim = 1;
-			NonLocMean->data.resize(Nx[0],Ny[0],Nz[0]);
+            NonLocMean->data.resize(Nx[i],Ny[i],Nz[i]);
-			meshData[0].vars.push_back(NonLocMean);
+            meshData[i].vars.push_back(NonLocMean);
-			fillDouble[0]->copy( NonLocalMean[0], NonLocMean->data );
+            fillDouble[i]->copy( NonLocalMean[i], NonLocMean->data );
-			std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
+            // Segmented Data
-			SegData->name = "Segmented Data";
+            auto SegData = std::make_shared<IO::Variable>();
-			SegData->type = IO::VariableType::VolumeVariable;
+            SegData->name = "Segmented_Data_" + std::to_string( i );
-			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 );
-
-			/*// Segmented 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[i],Ny[i],Nz[i]);
             meshData[i].vars.push_back(SegData);
             fillDouble[i]->copy( ID[i], SegData->data );
             // Signed Distance
-        std::shared_ptr<IO::Variable> DistData( new IO::Variable() );
+            auto DistData = std::make_shared<IO::Variable>();
-	    DistData->name = "Signed Distance";
+            DistData->name = "Signed_Distance_" + std::to_string( i );
             DistData->type = IO::VariableType::VolumeVariable;
             DistData->dim = 1;
             DistData->data.resize(Nx[i],Ny[i],Nz[i]);
             meshData[i].vars.push_back(DistData);
             fillDouble[i]->copy( Dist[i], DistData->data );
             // Smoothed Data
-        std::shared_ptr<IO::Variable> SmoothData( new IO::Variable() );
+            auto SmoothData = std::make_shared<IO::Variable>();
-	    SmoothData->name = "Smoothed Data";
+            SmoothData->name = "Smoothed_Data_" + std::to_string( i );
             SmoothData->type = IO::VariableType::VolumeVariable;
             SmoothData->dim = 1;
             SmoothData->data.resize(Nx[i],Ny[i],Nz[i]);
@@ -444,20 +398,18 @@ int main(int argc, char **argv)
             meshData[0].vars.push_back(filter_Dist2_var);
             fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data );
         #endif
-			 */
         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]));
+        auto Averages = std::make_shared<Minkowski>(Dm[0]);
         
-		Array <char> phase_label(Nx[0],Ny[0],Nz[0]);
+        Array <char> phase_label(Nx[0]+2,Ny[0]+2,Nz[0]+2);
-		Array <double> phase_distance(Nx[0],Ny[0],Nz[0]);
+        Array <double> phase_distance(Nx[0]+2,Ny[0]+2,Nz[0]+2);
         // Analyze the wetting fluid
         for (int k=1;k<Nz[0]+1;k++) {
             for (int j=1;j<Ny[0]+1;j++) {