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

common/Array.hpp

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

common/Communication.hpp

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

common/MPI_Helpers.hpp

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

common/Utilities.cpp

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

tests/lbpm_color_simulator.cpp

@@ -7,6 +7,7 @@
 #include <fstream>
 
 #include "common/ScaLBL.h"
+#include "common/Communication.h"
 #include "common/TwoPhase.h"
 #include "common/MPI_Helpers.h"
 
@@ -706,6 +707,37 @@ int main(int argc, char **argv)
         ThreadPool::setProcessAffinity(procs);
     }
     ThreadPool tpool(N_threads);
+
+    // Create the MeshDataStruct
+    fillHalo<double> fillData(Dm.Comm,Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
+    std::vector<IO::MeshDataStruct> meshData(1);
+    meshData[0].meshName = "domain";
+    meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
+    std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
+    std::shared_ptr<IO::Variable> PressVar( new IO::Variable() );
+    std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
+    std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
+    PhaseVar->name = "phase";
+    PhaseVar->type = IO::VolumeVariable;
+    PhaseVar->dim = 1;
+    PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(PhaseVar);
+    PressVar->name = "Pressure";
+    PressVar->type = IO::VolumeVariable;
+    PressVar->dim = 1;
+    PressVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(PressVar);
+    SignDistVar->name = "SignDist";
+    SignDistVar->type = IO::VolumeVariable;
+    SignDistVar->dim = 1;
+    SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(SignDistVar);
+    BlobIDVar->name = "BlobID";
+    BlobIDVar->type = IO::VolumeVariable;
+    BlobIDVar->dim = 1;
+    BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
+    meshData[0].vars.push_back(BlobIDVar);
+
 	//************ MAIN ITERATION LOOP ***************************************/
     PROFILE_START("Loop");
 	int timestep = -1;
@@ -714,6 +746,7 @@ int main(int argc, char **argv)
     writeIDMap(ID_map_struct(),0,id_map_filename);
     AnalysisWaitIdStruct work_ids;
 	while (timestep < timestepMax && err > tol ) {
+        if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
         PROFILE_START("Update");
 
 		//*************************************************************************
@@ -837,8 +870,8 @@ int main(int argc, char **argv)
         // Run the analysis, blob identification, and write restart files
         run_analysis(timestep,RESTART_INTERVAL,rank_info,*Averages,last_ids,last_index,last_id_map,
             Nx,Ny,Nz,pBC,beta,err,Phi,Pressure,Velocity,ID,f_even,f_odd,Den,
-            LocalRestartFile,tpool,work_ids);
+            LocalRestartFile,meshData,fillData,tpool,work_ids);
-
+        PROFILE_SAVE("lbpm_color_simulator",false);
 	}
     tpool.wait_pool_finished();
     PROFILE_STOP("Loop");

tests/lbpm_color_simulator.h

@@ -1,21 +1,26 @@
 // Run the analysis, blob identification, and write restart files
 #include "common/Array.h"
 #include "common/Communication.h"
+#include "common/MPI_Helpers.h"
+#include "IO/MeshDatabase.h"
 
+//#define ANALYSIS_INTERVAL 6
 #define ANALYSIS_INTERVAL 1000
 #define BLOBID_INTERVAL 250
 
 enum AnalysisType{ AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02, 
-    CopyAverages=0x04, CalcDist=0x08, CreateRestart=0x10 };
+    CopyAverages=0x04, CalcDist=0x08, CreateRestart=0x10, WriteVis=0x20 };
 
 
 // Structure used to store ids
 struct AnalysisWaitIdStruct {
     ThreadPool::thread_id_t blobID;
     ThreadPool::thread_id_t analysis;
+    ThreadPool::thread_id_t vis;
     ThreadPool::thread_id_t restart;
 };
 
+
 // Helper class to write the restart file from a seperate thread
 class WriteRestartWorkItem: public ThreadPool::WorkItem
 {
@@ -28,8 +33,6 @@ public:
         PROFILE_START("Save Checkpoint",1);
         WriteCheckpoint(filename,cDen.get(),cDistEven.get(),cDistOdd.get(),N);
         PROFILE_STOP("Save Checkpoint",1);
-
-        PROFILE_SAVE("lbpm_color_simulator",1);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     };
 private:
@@ -39,6 +42,7 @@ private:
     const int N;
 };
 
+
 // Helper class to compute the blob ids
 static const std::string id_map_filename = "lbpm_id_map.txt";
 typedef std::shared_ptr<std::pair<int,IntArray> > BlobIDstruct;
@@ -122,71 +126,42 @@ private:
     BlobIDList new_list;
 };
 
+
+// Helper class to write the vis file from a thread
 class WriteVisWorkItem: public ThreadPool::WorkItem
 {
 public:
-	WriteVisWorkItem(AnalysisType type_, TwoPhase& Averages_):
+    WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
-	type(type_), Averages(Averages_) { }
+        TwoPhase& Avgerages_, fillHalo<double>& fillData_ ):
-
+        timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_) {}
     virtual void run() {
         ThreadPool::WorkItem::d_state = 1;  // Change state to in progress
+        PROFILE_START("Save Vis",1);
+        ASSERT(visData[0].vars[0]->name=="phase");
+        ASSERT(visData[0].vars[1]->name=="Pressure");
+        ASSERT(visData[0].vars[2]->name=="SignDist");
+        ASSERT(visData[0].vars[3]->name=="BlobID");
+        Array<double>& PhaseData = visData[0].vars[0]->data;
+        Array<double>& PressData = visData[0].vars[1]->data;
+        Array<double>& SignData  = visData[0].vars[2]->data;
+        Array<double>& BlobData  = visData[0].vars[3]->data;
+        fillData.copy(Averages.SDn,PhaseData);
+        fillData.copy(Averages.Press,PressData);
+        fillData.copy(Averages.SDs,SignData);
+        fillData.copy(Averages.Label_NWP,BlobData);
         MPI_Comm newcomm;
         MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
-        // Write  VisIT files
+        IO::writeData( timestep, visData, 2, newcomm );
-        PROFILE_START("Save Vis",1);
-        // Create the MeshDataStruct
-        Nx=Averages.Dm.Nx;
-        Ny=Averages.Dm.Ny;
-        Nz=Averages.Dm.Nz;
-        Lx=Averages.Dm.Lx;
-        Ly=Averages.Dm.Ly;
-        Lz=Averages.Dm.Lz;
-        fillHalo<double> fillData(newcomm,Averages.Dm.rank_info,Nx-2,Ny-2,Nz-2,1,1,1,0,1);
-        std::vector<IO::MeshDataStruct> meshData(1);
-        meshData[0].meshName = "domain";
-        meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Averages.Dm.rank_info,Nx-2,Ny-2,Nz-2,Lx,Ly,Lz) );
-        std::shared_ptr<IO::Variable> PhaseVar( new IO::Variable() );
-        std::shared_ptr<IO::Variable> PressVar( new IO::Variable() );
-        std::shared_ptr<IO::Variable> SignDistVar( new IO::Variable() );
-        std::shared_ptr<IO::Variable> BlobIDVar( new IO::Variable() );
-        PhaseVar->name = "phase";
-        PhaseVar->type = IO::VolumeVariable;
-        PhaseVar->dim = 1;
-        PhaseVar->data.resize(Nx-2,Ny-2,Nz-2);
-        meshData[0].vars.push_back(PhaseVar);
-        PressVar->name = "Pressure";
-        PressVar->type = IO::VolumeVariable;
-        PressVar->dim = 1;
-        PressVar->data.resize(Nx-2,Ny-2,Nz-2);
-        meshData[0].vars.push_back(PressVar);
-        SignDistVar->name = "SignDist";
-        SignDistVar->type = IO::VolumeVariable;
-        SignDistVar->dim = 1;
-        SignDistVar->data.resize(Nx-2,Ny-2,Nz-2);
-        meshData[0].vars.push_back(SignDistVar);
-        BlobIDVar->name = "BlobID";
-        BlobIDVar->type = IO::VolumeVariable;
-        BlobIDVar->dim = 1;
-        BlobIDVar->data.resize(Nx-2,Ny-2,Nz-2);
-        meshData[0].vars.push_back(BlobIDVar);
-
-        fillData.copy(Averages.SDn,PhaseVar->data);
-        fillData.copy(Averages.SDs,SignDistVar->data);
-        fillData.copy(Averages.Label_NWP,BlobIDVar->data);
-        IO::writeData( 0, meshData, 2, newcomm );
-
         MPI_Comm_free(&newcomm);
         PROFILE_STOP("Save Vis",1);
-
-        PROFILE_SAVE("lbpm_color_simulator",1);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     };
 private:
     WriteVisWorkItem();
-    AnalysisType type;
+    int timestep;
+    std::vector<IO::MeshDataStruct>& visData;
     TwoPhase& Averages;
-    int Nx,Ny,Nz;
+    fillHalo<double>& fillData;
-    double Lx,Ly,Lz;
 };
 
 
@@ -225,6 +200,7 @@ public:
             Averages.PrintComponents(timestep);
             PROFILE_STOP("Compute dist",1);
         }
+        PROFILE_SAVE("lbpm_color_simulator",false);
         ThreadPool::WorkItem::d_state = 2;  // Change state to finished
     }
 private:
@@ -237,6 +213,8 @@ private:
     double beta;
 };
 
+
+
 // Function to start the analysis
 void run_analysis( int timestep, int restart_interval, 
     const RankInfoStruct& rank_info, TwoPhase& Averages,
@@ -244,7 +222,8 @@ void run_analysis( int timestep, int restart_interval,
     int Nx, int Ny, int Nz, bool pBC, double beta, double err,
     const double *Phi, double *Pressure, const double *Velocity, 
     const char *ID, const double *f_even, const double *f_odd, const double *Den, 
-    const char *LocalRestartFile, ThreadPool& tpool, AnalysisWaitIdStruct& wait )
+    const char *LocalRestartFile, std::vector<IO::MeshDataStruct>& visData, fillHalo<double>& fillData,
+    ThreadPool& tpool, AnalysisWaitIdStruct& wait )
 {
     int N = Nx*Ny*Nz;
 
@@ -266,7 +245,6 @@ void run_analysis( int timestep, int restart_interval,
             type = static_cast<AnalysisType>( type | IdentifyBlobs );
         }
     #endif
-
     if ( timestep%ANALYSIS_INTERVAL == 0 ) {
         // Copy the averages to the CPU (and identify blobs)
         type = static_cast<AnalysisType>( type | CopyAverages );
@@ -280,6 +258,12 @@ void run_analysis( int timestep, int restart_interval,
         // Write the restart file
         type = static_cast<AnalysisType>( type | CreateRestart );
     }
+    if (timestep%restart_interval == 0) {
+        // Write the visualization data
+        type = static_cast<AnalysisType>( type | WriteVis );
+        type = static_cast<AnalysisType>( type | CopyAverages );
+        type = static_cast<AnalysisType>( type | IdentifyBlobs );
+    }
     
     // Return if we are not doing anything
     if ( type == AnalyzeNone )
@@ -305,14 +289,22 @@ void run_analysis( int timestep, int restart_interval,
     }
     if ( (type&CopyAverages) != 0 ) {
         // Copy the members of Averages to the cpu (phase was copied above)
+        // Wait 
+        PROFILE_START("Copy-Wait",1);
+        tpool.wait(wait.analysis);
+        tpool.wait(wait.vis);   // Make sure we are done using analysis before modifying
+        PROFILE_STOP("Copy-Wait",1);
+        PROFILE_START("Copy-Pressure",1);
         ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
         memcpy(Averages.Phase.get(),phase->get(),N*sizeof(double));
         DeviceBarrier();
+        PROFILE_STOP("Copy-Pressure",1);
+        PROFILE_START("Copy-Averages",1);
         CopyToHost(Averages.Press.get(),Pressure,N*sizeof(double));
         CopyToHost(Averages.Vel_x.get(),&Velocity[0],N*sizeof(double));
         CopyToHost(Averages.Vel_y.get(),&Velocity[N],N*sizeof(double));
         CopyToHost(Averages.Vel_z.get(),&Velocity[2*N],N*sizeof(double));
-
+        PROFILE_STOP("Copy-Averages",1);
     }
     std::shared_ptr<double> cDen, cDistEven, cDistOdd;
     if ( (type&CreateRestart) != 0 ) {
@@ -349,6 +341,7 @@ void run_analysis( int timestep, int restart_interval,
             type,timestep,Averages,last_index,last_id_map,beta);
         work->add_dependency(wait.blobID);
         work->add_dependency(wait.analysis);
+        work->add_dependency(wait.vis);     // Make sure we are done using analysis before modifying
         wait.analysis = tpool.add_work(work);
     }
 
@@ -362,12 +355,27 @@ void run_analysis( int timestep, int restart_interval,
         } else {
             // Not clear yet
         }
+        // Wait for previous restart files to finish writing (not necessary, but helps to ensure memory usage is limited)
+        tpool.wait(wait.restart);
         // Write the restart file (using a seperate thread)
         WriteRestartWorkItem *work = new WriteRestartWorkItem(LocalRestartFile,cDen,cDistEven,cDistOdd,N);
         work->add_dependency(wait.restart);
         wait.restart = tpool.add_work(work);
     }
+
+    // Save the results for visualization
+    if ( (type&CreateRestart) != 0 ) {
+        // Wait for previous restart files to finish writing (not necessary, but helps to ensure memory usage is limited)
+        tpool.wait(wait.vis);
+        // Write the vis files
+        ThreadPool::WorkItem *work = new WriteVisWorkItem( timestep, visData, Averages, fillData );
+        work->add_dependency(wait.blobID);
+        work->add_dependency(wait.analysis);
+        work->add_dependency(wait.vis);
+        wait.vis = tpool.add_work(work);
+    }
     PROFILE_STOP("start_analysis");
 }
 
 
+