From db1fe9e327f421dabcc12b1fe0719aba5029f241 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 30 Jul 2018 16:51:37 -0400 Subject: [PATCH] working on DECL isosurface tools --- analysis/decl.cpp | 417 ++++++++++--------- analysis/decl.h | 8 +- analysis/runAnalysis.cpp | 860 +++++++++++++++++++-------------------- 3 files changed, 649 insertions(+), 636 deletions(-) diff --git a/analysis/decl.cpp b/analysis/decl.cpp index af5349e0..e63c8289 100644 --- a/analysis/decl.cpp +++ b/analysis/decl.cpp @@ -33,55 +33,64 @@ Point Vertex::coords(unsigned long int idx){ return P; } + +Halfedge::Halfedge(){ +} + +Halfedge::~Halfedge(){ + +} unsigned long int Halfedge::v1(unsigned long int edge){ - return HalfEdge(0,edge); + return data(0,edge); } unsigned long int Halfedge::v2(unsigned long int edge){ - return HalfEdge(1,edge); + return data(1,edge); } unsigned long int Halfedge::face(unsigned long int edge){ - return HalfEdge(2,edge); + return data(2,edge); } unsigned long int Halfedge::twin(unsigned long int edge){ - return HalfEdge(3,edge); + return data(3,edge); } unsigned long int Halfedge::prev(unsigned long int edge){ - return HalfEdge(4,edge); + return data(4,edge); } unsigned long int Halfedge::next(unsigned long int edge){ - return HalfEdge(5,edge); + return data(5,edge); +} + +DECL::DECL(){ +} + +DECL::~DECL(){ + } void DECL::LocalIsosurface(const DoubleArray A, double value, int i, int j, int k){ Point P,Q; Point PlaceHolder; - double temp; Point C0,C1,C2,C3,C4,C5,C6,C7; - + int TriangleCount; - int NewVertexCount; + int NewVertexCount; int CubeIndex; int nTris, nVert; - + Point VertexList[12]; Point NewVertexList[12]; int LocalRemap[12]; - + DTMutableList cellvertices = DTMutableList(20); IntArray Triangles = IntArray(3,20); // Values from array 'A' at the cube corners double CubeValues[8]; - - int Nx = A.size(0); - int Ny = A.size(1); - int Nz = A.size(2); - + // Points corresponding to cube corners C0.x = 0.0; C0.y = 0.0; C0.z = 0.0; C1.x = 1.0; C1.y = 0.0; C1.z = 0.0; @@ -92,188 +101,186 @@ void DECL::LocalIsosurface(const DoubleArray A, double value, int i, int j, int C6.x = 1.0; C6.y = 1.0; C6.z = 1.0; C7.x = 0.0; C7.y = 1.0; C7.z = 1.0; - CubeValues[0] = A(i,j,k) - value; - CubeValues[1] = A(i+1,j,k) - value; - CubeValues[2] = A(i+1,j+1,k) - value; - CubeValues[3] = A(i,j+1,k) - value; - CubeValues[4] = A(i,j,k+1) - value; - CubeValues[5] = A(i+1,j,k+1) - value; - CubeValues[6] = A(i+1,j+1,k+1) - value; - CubeValues[7] = A(i,j+1,k+1) -value; + CubeValues[0] = A(i,j,k) - value; + CubeValues[1] = A(i+1,j,k) - value; + CubeValues[2] = A(i+1,j+1,k) - value; + CubeValues[3] = A(i,j+1,k) - value; + CubeValues[4] = A(i,j,k+1) - value; + CubeValues[5] = A(i+1,j,k+1) - value; + CubeValues[6] = A(i+1,j+1,k+1) - value; + CubeValues[7] = A(i,j+1,k+1) -value; - //Determine the index into the edge table which - //tells us which vertices are inside of the surface - CubeIndex = 0; - if (CubeValues[0] < 0.0f) CubeIndex |= 1; - if (CubeValues[1] < 0.0f) CubeIndex |= 2; - if (CubeValues[2] < 0.0f) CubeIndex |= 4; - if (CubeValues[3] < 0.0f) CubeIndex |= 8; - if (CubeValues[4] < 0.0f) CubeIndex |= 16; - if (CubeValues[5] < 0.0f) CubeIndex |= 32; - if (CubeValues[6] < 0.0f) CubeIndex |= 64; - if (CubeValues[7] < 0.0f) CubeIndex |= 128; + //Determine the index into the edge table which + //tells us which vertices are inside of the surface + CubeIndex = 0; + if (CubeValues[0] < 0.0f) CubeIndex |= 1; + if (CubeValues[1] < 0.0f) CubeIndex |= 2; + if (CubeValues[2] < 0.0f) CubeIndex |= 4; + if (CubeValues[3] < 0.0f) CubeIndex |= 8; + if (CubeValues[4] < 0.0f) CubeIndex |= 16; + if (CubeValues[5] < 0.0f) CubeIndex |= 32; + if (CubeValues[6] < 0.0f) CubeIndex |= 64; + if (CubeValues[7] < 0.0f) CubeIndex |= 128; - //Find the vertices where the surface intersects the cube - if (edgeTable[CubeIndex] & 1){ - P = VertexInterp(C0,C1,CubeValues[0],CubeValues[1]); - VertexList[0] = P; - Q = C0; - } - if (edgeTable[CubeIndex] & 2){ - P = VertexInterp(C1,C2,CubeValues[1],CubeValues[2]); - VertexList[1] = P; - Q = C1; - } - if (edgeTable[CubeIndex] & 4){ - P = VertexInterp(C2,C3,CubeValues[2],CubeValues[3]); - VertexList[2] = P; - Q = C2; - } - if (edgeTable[CubeIndex] & 8){ - P = VertexInterp(C3,C0,CubeValues[3],CubeValues[0]); - VertexList[3] = P; - Q = C3; - } - if (edgeTable[CubeIndex] & 16){ - P = VertexInterp(C4,C5,CubeValues[4],CubeValues[5]); - VertexList[4] = P; - Q = C4; - } - if (edgeTable[CubeIndex] & 32){ - P = VertexInterp(C5,C6,CubeValues[5],CubeValues[6]); - VertexList[5] = P; - Q = C5; - } - if (edgeTable[CubeIndex] & 64){ - P = VertexInterp(C6,C7,CubeValues[6],CubeValues[7]); - VertexList[6] = P; - Q = C6; - } - if (edgeTable[CubeIndex] & 128){ - P = VertexInterp(C7,C4,CubeValues[7],CubeValues[4]); - VertexList[7] = P; - Q = C7; - } - if (edgeTable[CubeIndex] & 256){ - P = VertexInterp(C0,C4,CubeValues[0],CubeValues[4]); - VertexList[8] = P; - Q = C0; - } - if (edgeTable[CubeIndex] & 512){ - P = VertexInterp(C1,C5,CubeValues[1],CubeValues[5]); - VertexList[9] = P; - Q = C1; - } - if (edgeTable[CubeIndex] & 1024){ - P = VertexInterp(C2,C6,CubeValues[2],CubeValues[6]); - VertexList[10] = P; - Q = C2; - } - if (edgeTable[CubeIndex] & 2048){ - P = VertexInterp(C3,C7,CubeValues[3],CubeValues[7]); - VertexList[11] = P; - Q = C3; - } + //Find the vertices where the surface intersects the cube + if (edgeTable[CubeIndex] & 1){ + P = VertexInterp(C0,C1,CubeValues[0],CubeValues[1]); + VertexList[0] = P; + Q = C0; + } + if (edgeTable[CubeIndex] & 2){ + P = VertexInterp(C1,C2,CubeValues[1],CubeValues[2]); + VertexList[1] = P; + Q = C1; + } + if (edgeTable[CubeIndex] & 4){ + P = VertexInterp(C2,C3,CubeValues[2],CubeValues[3]); + VertexList[2] = P; + Q = C2; + } + if (edgeTable[CubeIndex] & 8){ + P = VertexInterp(C3,C0,CubeValues[3],CubeValues[0]); + VertexList[3] = P; + Q = C3; + } + if (edgeTable[CubeIndex] & 16){ + P = VertexInterp(C4,C5,CubeValues[4],CubeValues[5]); + VertexList[4] = P; + Q = C4; + } + if (edgeTable[CubeIndex] & 32){ + P = VertexInterp(C5,C6,CubeValues[5],CubeValues[6]); + VertexList[5] = P; + Q = C5; + } + if (edgeTable[CubeIndex] & 64){ + P = VertexInterp(C6,C7,CubeValues[6],CubeValues[7]); + VertexList[6] = P; + Q = C6; + } + if (edgeTable[CubeIndex] & 128){ + P = VertexInterp(C7,C4,CubeValues[7],CubeValues[4]); + VertexList[7] = P; + Q = C7; + } + if (edgeTable[CubeIndex] & 256){ + P = VertexInterp(C0,C4,CubeValues[0],CubeValues[4]); + VertexList[8] = P; + Q = C0; + } + if (edgeTable[CubeIndex] & 512){ + P = VertexInterp(C1,C5,CubeValues[1],CubeValues[5]); + VertexList[9] = P; + Q = C1; + } + if (edgeTable[CubeIndex] & 1024){ + P = VertexInterp(C2,C6,CubeValues[2],CubeValues[6]); + VertexList[10] = P; + Q = C2; + } + if (edgeTable[CubeIndex] & 2048){ + P = VertexInterp(C3,C7,CubeValues[3],CubeValues[7]); + VertexList[11] = P; + Q = C3; + } - NewVertexCount=0; - for (int idx=0;idx<12;idx++) - LocalRemap[idx] = -1; + NewVertexCount=0; + for (int idx=0;idx<12;idx++) + LocalRemap[idx] = -1; - for (int idx=0;triTable[CubeIndex][idx]!=-1;idx++) - { - if(LocalRemap[triTable[CubeIndex][idx]] == -1) - { - NewVertexList[NewVertexCount] = VertexList[triTable[CubeIndex][idx]]; - LocalRemap[triTable[CubeIndex][idx]] = NewVertexCount; - NewVertexCount++; - } - } + for (int idx=0;triTable[CubeIndex][idx]!=-1;idx++) + { + if(LocalRemap[triTable[CubeIndex][idx]] == -1) + { + NewVertexList[NewVertexCount] = VertexList[triTable[CubeIndex][idx]]; + LocalRemap[triTable[CubeIndex][idx]] = NewVertexCount; + NewVertexCount++; + } + } - for (int idx=0;idxV2 - HalfEdge(0,idx_edge) = V1; // first vertex - HalfEdge(1,idx_edge) = V2; // second vertex - HalfEdge(2,idx_edge) = idx; // triangle - HalfEdge(3,idx_edge) = -1; // twin - HalfEdge(4,idx_edge) = idx_edge+2; // previous edge - HalfEdge(5,idx_edge) = idx_edge+1; // next edge - idx_edge++; - // second edge: V2->V3 - HalfEdge(0,idx_edge) = V2; // first vertex - HalfEdge(1,idx_edge) = V3; // second vertex - HalfEdge(2,idx_edge) = idx; // triangle - HalfEdge(3,idx_edge) = -1; // twin - HalfEdge(4,idx_edge) = idx_edge-1; // previous edge - HalfEdge(5,idx_edge) = idx_edge+1; // next edge - idx_edge++; - // third edge: V3->V1 - HalfEdge(0,idx_edge) = V3; // first vertex - HalfEdge(1,idx_edge) = V1; // second vertex - HalfEdge(2,idx_edge) = idx; // triangle - HalfEdge(3,idx_edge) = -1; // twin - HalfEdge(4,idx_edge) = idx_edge-1; // previous edge - HalfEdge(5,idx_edge) = idx_edge-2; // next edge - idx_edge++; - } - int EdgeCount=idx_edge; - for (int idx=0; idxV2 + halfedge.data(0,idx_edge) = V1; // first vertex + halfedge.data(1,idx_edge) = V2; // second vertex + halfedge.data(2,idx_edge) = idx; // triangle + halfedge.data(3,idx_edge) = -1; // twin + halfedge.data(4,idx_edge) = idx_edge+2; // previous edge + halfedge.data(5,idx_edge) = idx_edge+1; // next edge + idx_edge++; + // second edge: V2->V3 + halfedge.data(0,idx_edge) = V2; // first vertex + halfedge.data(1,idx_edge) = V3; // second vertex + halfedge.data(2,idx_edge) = idx; // triangle + halfedge.data(3,idx_edge) = -1; // twin + halfedge.data(4,idx_edge) = idx_edge-1; // previous edge + halfedge.data(5,idx_edge) = idx_edge+1; // next edge + idx_edge++; + // third edge: V3->V1 + halfedge.data(0,idx_edge) = V3; // first vertex + halfedge.data(1,idx_edge) = V1; // second vertex + halfedge.data(2,idx_edge) = idx; // triangle + halfedge.data(3,idx_edge) = -1; // twin + halfedge.data(4,idx_edge) = idx_edge-1; // previous edge + halfedge.data(5,idx_edge) = idx_edge-2; // next edge + idx_edge++; + } + int EdgeCount=idx_edge; + for (int idx=0; idx cellvertices = DTMutableList(20); IntArray Triangles = IntArray(3,20); // Values from array 'A' at the cube corners double CubeValues[8]; - + int Nx = A.size(0); int Ny = A.size(1); int Nz = A.size(2); - + // Points corresponding to cube corners C0.x = 0.0; C0.y = 0.0; C0.z = 0.0; C1.x = 1.0; C1.y = 0.0; C1.z = 0.0; @@ -352,7 +365,7 @@ void Isosurface(DoubleArray &A, const double &v) C5.x = 1.0; C5.y = 0.0; C5.z = 1.0; C6.x = 1.0; C6.y = 1.0; C6.z = 1.0; C7.x = 0.0; C7.y = 1.0; C7.z = 1.0; - + for (int k=1; k data; private: - Array HalfEdge; }; // DECL @@ -45,12 +45,12 @@ public: unsigned long int face(); Vertex vertex; Halfedge halfedge; - void AddCube(); // need a function to add new faces based on marching cubes surface + void LocalIsosurface(const DoubleArray A, double value, int i, int j, int k); double origin(int edge); double EdgeAngle(int edge); Point TriNormal(int edge); - + private: unsigned long int *face_data; }; diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 51c9934d..54c86302 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -16,23 +16,23 @@ AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs) { - lhs = static_cast( - static_cast::type>(lhs) | - static_cast::type>(rhs) - ); - return lhs; + lhs = static_cast( + static_cast::type>(lhs) | + static_cast::type>(rhs) + ); + return lhs; } bool matches( AnalysisType x, AnalysisType y ) { - return static_cast::type>(x) & - static_cast::type>(y) != 0; + return static_cast::type>(x) & + static_cast::type>(y) != 0; } template void DeleteArray( const TYPE *p ) { - delete [] p; + delete [] p; } @@ -40,32 +40,32 @@ void DeleteArray( const TYPE *p ) class WriteRestartWorkItem: public ThreadPool::WorkItemRet { public: - WriteRestartWorkItem( const char* filename_, std::shared_ptr cphi_, std::shared_ptr cfq_, int N_ ): - filename(filename_), cphi(cphi_), cfq(cfq_), N(N_) {} - virtual void run() { - PROFILE_START("Save Checkpoint",1); - double value; - ofstream File(filename,ios::binary); - for (int n=0; n cphi_, std::shared_ptr cfq_, int N_ ): + filename(filename_), cphi(cphi_), cfq(cfq_), N(N_) {} + virtual void run() { + PROFILE_START("Save Checkpoint",1); + double value; + ofstream File(filename,ios::binary); + for (int n=0; n cfq,cphi; - // const DoubleArray& phase; - //const DoubleArray& dist; - const int N; + WriteRestartWorkItem(); + const char* filename; + std::shared_ptr cfq,cphi; + // const DoubleArray& phase; + //const DoubleArray& dist; + const int N; }; @@ -74,78 +74,78 @@ static const std::string id_map_filename = "lbpm_id_map.txt"; class BlobIdentificationWorkItem1: public ThreadPool::WorkItemRet { public: - BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, - std::shared_ptr phase_, const DoubleArray& dist_, - BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper&& comm_ ): - timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), - phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_)) - { - } - ~BlobIdentificationWorkItem1() { } - virtual void run() { - // Compute the global blob id and compare to the previous version - PROFILE_START("Identify blobs",1); - double vF = 0.0; - double vS = -1.0; // one voxel buffer region around solid - IntArray& ids = new_index->second; - new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,comm.comm); - PROFILE_STOP("Identify blobs",1); - } + BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, + std::shared_ptr phase_, const DoubleArray& dist_, + BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper&& comm_ ): + timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), + phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_)) +{ +} + ~BlobIdentificationWorkItem1() { } + virtual void run() { + // Compute the global blob id and compare to the previous version + PROFILE_START("Identify blobs",1); + double vF = 0.0; + double vS = -1.0; // one voxel buffer region around solid + IntArray& ids = new_index->second; + new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,comm.comm); + PROFILE_STOP("Identify blobs",1); + } private: - BlobIdentificationWorkItem1(); - int timestep; - int Nx, Ny, Nz; - const RankInfoStruct& rank_info; - std::shared_ptr phase; - const DoubleArray& dist; - BlobIDstruct last_id, new_index, new_id; - BlobIDList new_list; - runAnalysis::commWrapper comm; + BlobIdentificationWorkItem1(); + int timestep; + int Nx, Ny, Nz; + const RankInfoStruct& rank_info; + std::shared_ptr phase; + const DoubleArray& dist; + BlobIDstruct last_id, new_index, new_id; + BlobIDList new_list; + runAnalysis::commWrapper comm; }; class BlobIdentificationWorkItem2: public ThreadPool::WorkItemRet { public: - BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, - std::shared_ptr phase_, const DoubleArray& dist_, - BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ , runAnalysis::commWrapper&& comm_ ): - timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), - phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_)) - { - } - ~BlobIdentificationWorkItem2() { } - virtual void run() { - // Compute the global blob id and compare to the previous version - PROFILE_START("Identify blobs maps",1); - const IntArray& ids = new_index->second; - static int max_id = -1; - new_id->first = new_index->first; - new_id->second = new_index->second; - if ( last_id.get()!=NULL ) { - // Compute the timestep-timestep map - const IntArray& old_ids = last_id->second; - ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,comm.comm); - // Renumber the current timestep's ids - getNewIDs(map,max_id,*new_list); - renumberIDs(*new_list,new_id->second); - writeIDMap(map,timestep,id_map_filename); - } else { - max_id = -1; - ID_map_struct map(new_id->first); - getNewIDs(map,max_id,*new_list); - writeIDMap(map,timestep,id_map_filename); - } - PROFILE_STOP("Identify blobs maps",1); - } + BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, + std::shared_ptr phase_, const DoubleArray& dist_, + BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ , runAnalysis::commWrapper&& comm_ ): + timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), + phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_)) +{ +} + ~BlobIdentificationWorkItem2() { } + virtual void run() { + // Compute the global blob id and compare to the previous version + PROFILE_START("Identify blobs maps",1); + const IntArray& ids = new_index->second; + static int max_id = -1; + new_id->first = new_index->first; + new_id->second = new_index->second; + if ( last_id.get()!=NULL ) { + // Compute the timestep-timestep map + const IntArray& old_ids = last_id->second; + ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,comm.comm); + // Renumber the current timestep's ids + getNewIDs(map,max_id,*new_list); + renumberIDs(*new_list,new_id->second); + writeIDMap(map,timestep,id_map_filename); + } else { + max_id = -1; + ID_map_struct map(new_id->first); + getNewIDs(map,max_id,*new_list); + writeIDMap(map,timestep,id_map_filename); + } + PROFILE_STOP("Identify blobs maps",1); + } private: - BlobIdentificationWorkItem2(); - int timestep; - int Nx, Ny, Nz; - const RankInfoStruct& rank_info; - std::shared_ptr phase; - const DoubleArray& dist; - BlobIDstruct last_id, new_index, new_id; - BlobIDList new_list; - runAnalysis::commWrapper comm; + BlobIdentificationWorkItem2(); + int timestep; + int Nx, Ny, Nz; + const RankInfoStruct& rank_info; + std::shared_ptr phase; + const DoubleArray& dist; + BlobIDstruct last_id, new_index, new_id; + BlobIDList new_list; + runAnalysis::commWrapper comm; }; @@ -153,36 +153,36 @@ private: class WriteVisWorkItem: public ThreadPool::WorkItemRet { public: - WriteVisWorkItem( int timestep_, std::vector& visData_, - TwoPhase& Avgerages_, fillHalo& fillData_, runAnalysis::commWrapper&& comm_ ): - timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_), comm(std::move(comm_)) - { - } - ~WriteVisWorkItem() { } - virtual void run() { - 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& PhaseData = visData[0].vars[0]->data; - Array& PressData = visData[0].vars[1]->data; - Array& SignData = visData[0].vars[2]->data; - Array& 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); - IO::writeData( timestep, visData, comm.comm ); - PROFILE_STOP("Save Vis",1); - }; + WriteVisWorkItem( int timestep_, std::vector& visData_, + TwoPhase& Avgerages_, fillHalo& fillData_, runAnalysis::commWrapper&& comm_ ): + timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_), comm(std::move(comm_)) +{ +} + ~WriteVisWorkItem() { } + virtual void run() { + 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& PhaseData = visData[0].vars[0]->data; + Array& PressData = visData[0].vars[1]->data; + Array& SignData = visData[0].vars[2]->data; + Array& 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); + IO::writeData( timestep, visData, comm.comm ); + PROFILE_STOP("Save Vis",1); + }; private: - WriteVisWorkItem(); - int timestep; - std::vector& visData; - TwoPhase& Averages; - fillHalo& fillData; - runAnalysis::commWrapper comm; + WriteVisWorkItem(); + int timestep; + std::vector& visData; + TwoPhase& Averages; + fillHalo& fillData; + runAnalysis::commWrapper comm; }; @@ -191,46 +191,46 @@ private: class AnalysisWorkItem: public ThreadPool::WorkItemRet { public: - AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, - BlobIDstruct ids, BlobIDList id_list_, double beta_ ): - type(type_), timestep(timestep_), Averages(Averages_), - blob_ids(ids), id_list(id_list_), beta(beta_) { } - ~AnalysisWorkItem() { } - virtual void run() { - Averages.NumberComponents_NWP = blob_ids->first; - Averages.Label_NWP = blob_ids->second; - Averages.Label_NWP_map = *id_list; - Averages.NumberComponents_WP = 1; - Averages.Label_WP.fill(0.0); - if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { - // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus); - } - if ( matches(type,AnalysisType::ComputeAverages) ) { - PROFILE_START("Compute dist",1); - Averages.Initialize(); - Averages.ComputeDelPhi(); - Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn); - Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus); - Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus); - Averages.UpdateMeshValues(); - Averages.ComputeLocal(); - Averages.Reduce(); - Averages.PrintAll(timestep); - Averages.Initialize(); - Averages.ComponentAverages(); - Averages.SortBlobs(); - Averages.PrintComponents(timestep); - PROFILE_STOP("Compute dist",1); - } - } + AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, + BlobIDstruct ids, BlobIDList id_list_, double beta_ ): + type(type_), timestep(timestep_), Averages(Averages_), + blob_ids(ids), id_list(id_list_), beta(beta_) { } + ~AnalysisWorkItem() { } + virtual void run() { + Averages.NumberComponents_NWP = blob_ids->first; + Averages.Label_NWP = blob_ids->second; + Averages.Label_NWP_map = *id_list; + Averages.NumberComponents_WP = 1; + Averages.Label_WP.fill(0.0); + if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { + // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus); + } + if ( matches(type,AnalysisType::ComputeAverages) ) { + PROFILE_START("Compute dist",1); + Averages.Initialize(); + Averages.ComputeDelPhi(); + Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn); + Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus); + Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus); + Averages.UpdateMeshValues(); + Averages.ComputeLocal(); + Averages.Reduce(); + Averages.PrintAll(timestep); + Averages.Initialize(); + Averages.ComponentAverages(); + Averages.SortBlobs(); + Averages.PrintComponents(timestep); + PROFILE_STOP("Compute dist",1); + } + } private: - AnalysisWorkItem(); - AnalysisType type; - int timestep; - TwoPhase& Averages; - BlobIDstruct blob_ids; - BlobIDList id_list; - double beta; + AnalysisWorkItem(); + AnalysisType type; + int timestep; + TwoPhase& Averages; + BlobIDstruct blob_ids; + BlobIDList id_list; + double beta; }; @@ -238,44 +238,44 @@ private: * MPI comm wrapper for use with analysis * ******************************************************************/ runAnalysis::commWrapper::commWrapper( int tag_, MPI_Comm comm_, runAnalysis* analysis_ ): - comm(comm_), - tag(tag_), - analysis(analysis_) + comm(comm_), + tag(tag_), + analysis(analysis_) { } runAnalysis::commWrapper::commWrapper( commWrapper &&rhs ): - comm(rhs.comm), - tag(rhs.tag), - analysis(rhs.analysis) + comm(rhs.comm), + tag(rhs.tag), + analysis(rhs.analysis) { - rhs.tag = -1; + rhs.tag = -1; } runAnalysis::commWrapper::~commWrapper() { - if ( tag == -1 ) - return; - MPI_Barrier( comm ); - analysis->d_comm_used[tag] = false; + if ( tag == -1 ) + return; + MPI_Barrier( comm ); + analysis->d_comm_used[tag] = false; } runAnalysis::commWrapper runAnalysis::getComm( ) { - // Get a tag from root - int tag = -1; - if ( d_rank == 0 ) { - for (int i=0; i<1024; i++) { - if ( !d_comm_used[i] ) { - tag = i; - break; - } - } - if ( tag == -1 ) - ERROR("Unable to get comm"); - } - MPI_Bcast( &tag, 1, MPI_INT, 0, d_comm ); - d_comm_used[tag] = true; - if ( d_comms[tag] == MPI_COMM_NULL ) - MPI_Comm_dup( MPI_COMM_WORLD, &d_comms[tag] ); - return commWrapper(tag,d_comms[tag],this); + // Get a tag from root + int tag = -1; + if ( d_rank == 0 ) { + for (int i=0; i<1024; i++) { + if ( !d_comm_used[i] ) { + tag = i; + break; + } + } + if ( tag == -1 ) + ERROR("Unable to get comm"); + } + MPI_Bcast( &tag, 1, MPI_INT, 0, d_comm ); + d_comm_used[tag] = true; + if ( d_comms[tag] == MPI_COMM_NULL ) + MPI_Comm_dup( MPI_COMM_WORLD, &d_comms[tag] ); + return commWrapper(tag,d_comms[tag],this); } @@ -283,29 +283,29 @@ runAnalysis::commWrapper runAnalysis::getComm( ) * Constructor/Destructors * ******************************************************************/ runAnalysis::runAnalysis( std::shared_ptr db, - const RankInfoStruct& rank_info, std::shared_ptr ScaLBL_Comm, std::shared_ptr Dm, - int Np, bool Regular, double beta, IntArray Map ): - d_Np( Np ), - d_beta( beta ), - 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) + const RankInfoStruct& rank_info, std::shared_ptr ScaLBL_Comm, std::shared_ptr Dm, + int Np, bool Regular, double beta, IntArray Map ): + d_Np( Np ), + d_beta( beta ), + 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) { char rankString[20]; sprintf(rankString,"%05d",Dm->rank()); - d_N[0] = Dm->Nx; - d_N[1] = Dm->Ny; - d_N[2] = Dm->Nz; - d_restart_interval = db->getScalar( "restart_interval" ); - d_analysis_interval = db->getScalar( "analysis_interval" ); - d_blobid_interval = db->getScalar( "blobid_interval" ); - d_visualization_interval = db->getScalar( "visualization_interval" ); - auto restart_file = db->getScalar( "restart_file" ); - d_restartFile = restart_file + "." + rankString; - d_rank = MPI_WORLD_RANK(); + d_N[0] = Dm->Nx; + d_N[1] = Dm->Ny; + d_N[2] = Dm->Nz; + d_restart_interval = db->getScalar( "restart_interval" ); + d_analysis_interval = db->getScalar( "analysis_interval" ); + d_blobid_interval = db->getScalar( "blobid_interval" ); + d_visualization_interval = db->getScalar( "visualization_interval" ); + auto restart_file = db->getScalar( "restart_file" ); + d_restartFile = restart_file + "." + rankString; + d_rank = MPI_WORLD_RANK(); writeIDMap(ID_map_struct(),0,id_map_filename); // Initialize IO for silo IO::initialize("","silo","false"); @@ -337,41 +337,41 @@ runAnalysis::runAnalysis( std::shared_ptr db, BlobIDVar->dim = 1; BlobIDVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); d_meshData[0].vars.push_back(BlobIDVar); - // Initialize the comms - MPI_Comm_dup(MPI_COMM_WORLD,&d_comm); - for (int i=0; i<1024; i++) { - d_comms[i] = MPI_COMM_NULL; - d_comm_used[i] = false; - } - // Initialize the threads - int N_threads = db->getWithDefault( "N_threads", 4 ); - auto method = db->getWithDefault( "load_balance", "default" ); - createThreads( method, N_threads ); + // Initialize the comms + MPI_Comm_dup(MPI_COMM_WORLD,&d_comm); + for (int i=0; i<1024; i++) { + d_comms[i] = MPI_COMM_NULL; + d_comm_used[i] = false; + } + // Initialize the threads + int N_threads = db->getWithDefault( "N_threads", 4 ); + auto method = db->getWithDefault( "load_balance", "default" ); + createThreads( method, N_threads ); } runAnalysis::~runAnalysis( ) { - // Finish processing analysis - finish(); - // Clear internal data - MPI_Comm_free( &d_comm ); - for (int i=0; i<1024; i++) { - if ( d_comms[i] != MPI_COMM_NULL ) - MPI_Comm_free(&d_comms[i]); - } + // Finish processing analysis + finish(); + // Clear internal data + MPI_Comm_free( &d_comm ); + for (int i=0; i<1024; i++) { + if ( d_comms[i] != MPI_COMM_NULL ) + MPI_Comm_free(&d_comms[i]); + } } void runAnalysis::finish( ) { - PROFILE_START("finish"); - // Wait for the work items to finish + PROFILE_START("finish"); + // Wait for the work items to finish d_tpool.wait_pool_finished(); - // Clear the wait ids - d_wait_blobID.reset(); - d_wait_analysis.reset(); - d_wait_vis.reset(); - d_wait_restart.reset(); - // Syncronize - MPI_Barrier( d_comm ); - PROFILE_STOP("finish"); + // Clear the wait ids + d_wait_blobID.reset(); + d_wait_analysis.reset(); + d_wait_vis.reset(); + d_wait_restart.reset(); + // Syncronize + MPI_Barrier( d_comm ); + PROFILE_STOP("finish"); } @@ -380,50 +380,50 @@ void runAnalysis::finish( ) ******************************************************************/ void print( const std::vector& ids ) { - if ( ids.empty() ) - return; - printf("%i",ids[0]); - for (size_t i=1; i0 && timestep%50==0 ) { // Keep a few blob identifications queued up to keep the processors busy, // allowing us to track the blobs as fast as possible @@ -450,28 +450,28 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep ) type |= AnalysisType::IdentifyBlobs; } #endif */ - if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) { - // Copy the averages to the CPU (and identify blobs) - //printf("Copy sim state, timestep=%i \n",timestep); - type |= AnalysisType::CopySimState; - type |= AnalysisType::IdentifyBlobs; - } - if ( timestep%d_analysis_interval == 0 ) { - // Run the analysis - //printf("Compute averages, timestep=%i \n",timestep); - type |= AnalysisType::ComputeAverages; - } - if (timestep%d_restart_interval == 0) { - // Write the restart file - type |= AnalysisType::CreateRestart; - } - if (timestep%d_visualization_interval == 0) { - // Write the visualization data - type |= AnalysisType::WriteVis; - type |= AnalysisType::CopySimState; - type |= AnalysisType::IdentifyBlobs; - } - return type; + if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) { + // Copy the averages to the CPU (and identify blobs) + //printf("Copy sim state, timestep=%i \n",timestep); + type |= AnalysisType::CopySimState; + type |= AnalysisType::IdentifyBlobs; + } + if ( timestep%d_analysis_interval == 0 ) { + // Run the analysis + //printf("Compute averages, timestep=%i \n",timestep); + type |= AnalysisType::ComputeAverages; + } + if (timestep%d_restart_interval == 0) { + // Write the restart file + type |= AnalysisType::CreateRestart; + } + if (timestep%d_visualization_interval == 0) { + // Write the visualization data + type |= AnalysisType::WriteVis; + type |= AnalysisType::CopySimState; + type |= AnalysisType::IdentifyBlobs; + } + return type; } @@ -480,28 +480,28 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep ) * Run the analysis * ******************************************************************/ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi, - double *Pressure, double *Velocity, double *fq, double *Den) + double *Pressure, double *Velocity, double *fq, double *Den) { - int N = d_N[0]*d_N[1]*d_N[2]; + int N = d_N[0]*d_N[1]*d_N[2]; - // Check which analysis steps we need to perform - auto type = computeAnalysisType( timestep ); - if ( type == AnalysisType::AnalyzeNone ) - return; + // Check which analysis steps we need to perform + auto type = computeAnalysisType( timestep ); + if ( type == AnalysisType::AnalyzeNone ) + return; - // Check how may queued items we have - if ( d_tpool.N_queued() > 20 ) { - std::cerr << "Analysis queue is getting behind, waiting ...\n"; - finish(); - } + // Check how may queued items we have + if ( d_tpool.N_queued() > 20 ) { + std::cerr << "Analysis queue is getting behind, waiting ...\n"; + finish(); + } - PROFILE_START("run"); + PROFILE_START("run"); - // Copy the appropriate variables to the host (so we can spawn new threads) - ScaLBL_DeviceBarrier(); - PROFILE_START("Copy data to host",1); - std::shared_ptr phase; - /* if ( matches(type,AnalysisType::CopyPhaseIndicator) || + // Copy the appropriate variables to the host (so we can spawn new threads) + ScaLBL_DeviceBarrier(); + PROFILE_START("Copy data to host",1); + std::shared_ptr phase; + /* if ( matches(type,AnalysisType::CopyPhaseIndicator) || matches(type,AnalysisType::ComputeAverages) || matches(type,AnalysisType::CopySimState) || matches(type,AnalysisType::IdentifyBlobs) ) @@ -531,113 +531,113 @@ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi, } delete [] TmpDat; } - */ - //if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { - if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) { - if (d_regular) - d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tplus); - else - ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double)); - //memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double)); - } - if ( timestep%d_analysis_interval == 0 ) { - if (d_regular) - d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tminus); - else - ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double)); - //memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double)); - } - //if ( matches(type,AnalysisType::CopySimState) ) { - if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) { - // Copy the members of Averages to the cpu (phase was copied above) - PROFILE_START("Copy-Pressure",1); - ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np); - ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np); - ScaLBL_DeviceBarrier(); - PROFILE_STOP("Copy-Pressure",1); - PROFILE_START("Copy-Wait",1); - PROFILE_STOP("Copy-Wait",1); - PROFILE_START("Copy-State",1); - //memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double)); - if (d_regular) - d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase); - else - ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double)); - // copy other variables - d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press); - d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x); - d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y); - d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z); - PROFILE_STOP("Copy-State",1); - } - std::shared_ptr cfq,cPhi; - //if ( matches(type,AnalysisType::CreateRestart) ) { - if (timestep%d_restart_interval==0){ - // Copy restart data to the CPU - cPhi = std::shared_ptr(new double[d_Np],DeleteArray); - cfq = std::shared_ptr(new double[19*d_Np],DeleteArray); - ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double)); - ScaLBL_CopyToHost(cPhi.get(),Phi,d_Np*sizeof(double)); - } - PROFILE_STOP("Copy data to host",1); + */ + //if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { + if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) { + if (d_regular) + d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tplus); + else + ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double)); + //memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double)); + } + if ( timestep%d_analysis_interval == 0 ) { + if (d_regular) + d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tminus); + else + ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double)); + //memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double)); + } + //if ( matches(type,AnalysisType::CopySimState) ) { + if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) { + // Copy the members of Averages to the cpu (phase was copied above) + PROFILE_START("Copy-Pressure",1); + ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np); + ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np); + ScaLBL_DeviceBarrier(); + PROFILE_STOP("Copy-Pressure",1); + PROFILE_START("Copy-Wait",1); + PROFILE_STOP("Copy-Wait",1); + PROFILE_START("Copy-State",1); + //memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double)); + if (d_regular) + d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase); + else + ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double)); + // copy other variables + d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press); + d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x); + d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y); + d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z); + PROFILE_STOP("Copy-State",1); + } + std::shared_ptr cfq,cPhi; + //if ( matches(type,AnalysisType::CreateRestart) ) { + if (timestep%d_restart_interval==0){ + // Copy restart data to the CPU + cPhi = std::shared_ptr(new double[d_Np],DeleteArray); + cfq = std::shared_ptr(new double[19*d_Np],DeleteArray); + ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double)); + ScaLBL_CopyToHost(cPhi.get(),Phi,d_Np*sizeof(double)); + } + PROFILE_STOP("Copy data to host",1); - // Spawn threads to do blob identification work - if ( matches(type,AnalysisType::IdentifyBlobs) ) { - phase = std::shared_ptr(new DoubleArray(d_N[0],d_N[1],d_N[2])); - d_ScaLBL_Comm->RegularLayout(d_Map,Phi,*phase); + // Spawn threads to do blob identification work + if ( matches(type,AnalysisType::IdentifyBlobs) ) { + phase = std::shared_ptr(new DoubleArray(d_N[0],d_N[1],d_N[2])); + d_ScaLBL_Comm->RegularLayout(d_Map,Phi,*phase); - BlobIDstruct new_index(new std::pair(0,IntArray())); - BlobIDstruct new_ids(new std::pair(0,IntArray())); - BlobIDList new_list(new std::vector()); - auto work1 = new BlobIdentificationWorkItem1(timestep,d_N[0],d_N[1],d_N[2],d_rank_info, - phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm()); - auto work2 = new BlobIdentificationWorkItem2(timestep,d_N[0],d_N[1],d_N[2],d_rank_info, - phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm()); - work1->add_dependency(d_wait_blobID); - work2->add_dependency(d_tpool.add_work(work1)); - d_wait_blobID = d_tpool.add_work(work2); - d_last_index = new_index; - d_last_ids = new_ids; - d_last_id_map = new_list; - } + BlobIDstruct new_index(new std::pair(0,IntArray())); + BlobIDstruct new_ids(new std::pair(0,IntArray())); + BlobIDList new_list(new std::vector()); + auto work1 = new BlobIdentificationWorkItem1(timestep,d_N[0],d_N[1],d_N[2],d_rank_info, + phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm()); + auto work2 = new BlobIdentificationWorkItem2(timestep,d_N[0],d_N[1],d_N[2],d_rank_info, + phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm()); + work1->add_dependency(d_wait_blobID); + work2->add_dependency(d_tpool.add_work(work1)); + d_wait_blobID = d_tpool.add_work(work2); + d_last_index = new_index; + d_last_ids = new_ids; + d_last_id_map = new_list; + } - // Spawn threads to do the analysis work - //if (timestep%d_restart_interval==0){ - // if ( matches(type,AnalysisType::ComputeAverages) ) { - if ( timestep%d_analysis_interval == 0 ) { - auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta); - work->add_dependency(d_wait_blobID); - work->add_dependency(d_wait_analysis); - work->add_dependency(d_wait_vis); // Make sure we are done using analysis before modifying - d_wait_analysis = d_tpool.add_work(work); - } + // Spawn threads to do the analysis work + //if (timestep%d_restart_interval==0){ + // if ( matches(type,AnalysisType::ComputeAverages) ) { + if ( timestep%d_analysis_interval == 0 ) { + auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta); + work->add_dependency(d_wait_blobID); + work->add_dependency(d_wait_analysis); + work->add_dependency(d_wait_vis); // Make sure we are done using analysis before modifying + d_wait_analysis = d_tpool.add_work(work); + } - // Spawn a thread to write the restart file - // if ( matches(type,AnalysisType::CreateRestart) ) { - if (timestep%d_restart_interval==0){ + // Spawn a thread to write the restart file + // if ( matches(type,AnalysisType::CreateRestart) ) { + if (timestep%d_restart_interval==0){ - if (d_rank==0) { - FILE *Rst = fopen("Restart.txt","w"); - fprintf(Rst,"%i\n",timestep+4); - fclose(Rst); - } - // Write the restart file (using a seperate thread) - auto work = new WriteRestartWorkItem(d_restartFile.c_str(),cPhi,cfq,d_Np); - work->add_dependency(d_wait_restart); - d_wait_restart = d_tpool.add_work(work); - } + if (d_rank==0) { + FILE *Rst = fopen("Restart.txt","w"); + fprintf(Rst,"%i\n",timestep+4); + fclose(Rst); + } + // Write the restart file (using a seperate thread) + auto work = new WriteRestartWorkItem(d_restartFile.c_str(),cPhi,cfq,d_Np); + work->add_dependency(d_wait_restart); + d_wait_restart = d_tpool.add_work(work); + } - // Save the results for visualization - // if ( matches(type,AnalysisType::CreateRestart) ) { - if (timestep%d_restart_interval==0){ - // Write the vis files - auto work = new WriteVisWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() ); - work->add_dependency(d_wait_blobID); - work->add_dependency(d_wait_analysis); - work->add_dependency(d_wait_vis); - d_wait_vis = d_tpool.add_work(work); - } - PROFILE_STOP("run"); + // Save the results for visualization + // if ( matches(type,AnalysisType::CreateRestart) ) { + if (timestep%d_restart_interval==0){ + // Write the vis files + auto work = new WriteVisWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() ); + work->add_dependency(d_wait_blobID); + work->add_dependency(d_wait_analysis); + work->add_dependency(d_wait_vis); + d_wait_vis = d_tpool.add_work(work); + } + PROFILE_STOP("run"); }