diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 7b466a24..f2597447 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -28,8 +28,15 @@ SubPhase::SubPhase(std::shared_ptr dm): //......................................... if (Dm->rank()==0){ + bool WriteHeader=false; + SUBPHASE = fopen("subphase.csv","r"); + if (SUBPHASE != NULL) + fclose(SUBPHASE); + else + WriteHeader=true; + SUBPHASE = fopen("subphase.csv","a+"); - if (ftell(SUBPHASE) == 0) + if (WriteHeader) { // If timelog is empty, write a short header to list the averages //fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n"); @@ -141,11 +148,12 @@ void SubPhase::Basic(){ if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4; imin=jmin=1; - // If inlet layers exist use these as default + // If inlet/outlet layers exist use these as default if (Dm->inlet_layers_x > 0) imin = Dm->inlet_layers_x; if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y; if (Dm->inlet_layers_z > 0) kmin = Dm->inlet_layers_z; - + if (Dm->outlet_layers_z > 0) kmax = Dm->outlet_layers_z; + nb.reset(); wb.reset(); /* //Dm->CommunicateMeshHalo(Phi); diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 06f83171..62802aa6 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -389,9 +389,9 @@ public: // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus); } if ( matches(type,AnalysisType::ComputeAverages) ) { - PROFILE_START("Compute subphase",1); + PROFILE_START("Compute basic averages",1); Averages.Basic(); - PROFILE_STOP("Compute subphase",1); + PROFILE_STOP("Compute basic averages",1); } } private: @@ -402,6 +402,29 @@ private: double beta; }; +class SubphaseWorkItem: public ThreadPool::WorkItemRet +{ +public: + SubphaseWorkItem( AnalysisType type_, int timestep_, SubPhase& Averages_ ): + type(type_), timestep(timestep_), Averages(Averages_){ } + ~SubphaseWorkItem() { } + virtual void run() { + + if ( matches(type,AnalysisType::ComputeAverages) ) { + PROFILE_START("Compute subphase",1); + Averages.Full(); + PROFILE_STOP("Compute subphase",1); + } + } +private: + SubphaseWorkItem(); + AnalysisType type; + int timestep; + SubPhase& Averages; + double beta; +}; + + /****************************************************************** * MPI comm wrapper for use with analysis * @@ -468,6 +491,7 @@ runAnalysis::runAnalysis( std::shared_ptr db, ThreadPool::thread_id_t d_wait_analysis; ThreadPool::thread_id_t d_wait_vis; ThreadPool::thread_id_t d_wait_restart; + ThreadPool::thread_id_t d_wait_subphase; char rankString[20]; sprintf(rankString,"%05d",Dm->rank()); @@ -479,6 +503,12 @@ runAnalysis::runAnalysis( std::shared_ptr db, d_blobid_interval = db->getScalar( "blobid_interval" ); d_visualization_interval = db->getScalar( "visualization_interval" ); auto restart_file = db->getScalar( "restart_file" ); + if (db->keyExists( "subphase_analysis_interval" )){ + d_subphase_analysis_interval = db->getScalar( "subphase_analysis_interval" ); + } + else{ + d_subphase_analysis_interval = INT_MAX; + } d_restartFile = restart_file + "." + rankString; d_rank = MPI_WORLD_RANK(); writeIDMap(ID_map_struct(),0,id_map_filename); @@ -892,9 +922,16 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do // if ( matches(type,AnalysisType::ComputeAverages) ) { if ( timestep%d_analysis_interval == 0 ) { auto work = new BasicWorkItem(type,timestep,Averages); - work->add_dependency(d_wait_analysis); // Make sure we are done using analysis before modifying + work->add_dependency(d_wait_subphase); // Make sure we are done using analysis before modifying + work->add_dependency(d_wait_analysis); d_wait_analysis = d_tpool.add_work(work); } + + if ( timestep%d_subphase_analysis_interval == 0 ) { + auto work = new BasicWorkItem(type,timestep,Averages); + work->add_dependency(d_wait_subphase); // Make sure we are done using analysis before modifying + d_wait_subphase = d_tpool.add_work(work); + } if (timestep%d_restart_interval==0){ std::shared_ptr cfq,cDen; diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h index de2c1053..d9c96e6b 100644 --- a/analysis/runAnalysis.h +++ b/analysis/runAnalysis.h @@ -7,7 +7,7 @@ #include "common/Communication.h" #include "common/ScaLBL.h" #include "threadpool/thread_pool.h" - +#include typedef std::shared_ptr> BlobIDstruct; typedef std::shared_ptr> BlobIDList; @@ -88,6 +88,7 @@ private: int d_Np; int d_rank; int d_restart_interval, d_analysis_interval, d_blobid_interval, d_visualization_interval; + int d_subphase_analysis_interval; double d_beta; bool d_regular; ThreadPool d_tpool; @@ -107,6 +108,7 @@ private: // 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_subphase; ThreadPool::thread_id_t d_wait_vis; ThreadPool::thread_id_t d_wait_restart; diff --git a/common/Domain.cpp b/common/Domain.cpp index 27a2359a..e88a6b25 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -77,6 +77,7 @@ Domain::Domain( std::shared_ptr db, MPI_Comm Communicator): Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0), Comm(MPI_COMM_NULL), inlet_layers_x(0), inlet_layers_y(0), inlet_layers_z(0), + outlet_layers_x(0), outlet_layers_y(0), outlet_layers_z(0), sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0), sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0), sendCount_xY(0), sendCount_yZ(0), sendCount_Xz(0), sendCount_XY(0), sendCount_YZ(0), sendCount_XZ(0), @@ -126,6 +127,12 @@ void Domain::initialize( std::shared_ptr db ) inlet_layers_y = InletCount[1]; inlet_layers_z = InletCount[2]; } + if (d_db->keyExists( "OutletLayers" )){ + auto OutletCount = d_db->getVector( "OutletLayers" ); + outlet_layers_x = OutletCount[0]; + outlet_layers_y = OutletCount[1]; + outlet_layers_z = OutletCount[2]; + } ASSERT( n.size() == 3u ); ASSERT( L.size() == 3u ); @@ -144,9 +151,13 @@ void Domain::initialize( std::shared_ptr db ) MPI_Comm_rank( Comm, &myrank ); rank_info = RankInfoStruct(myrank,nproc[0],nproc[1],nproc[2]); // inlet layers only apply to lower part of domain - if (nproc[0] > 0) inlet_layers_x = 0; - if (nproc[1] > 0) inlet_layers_y = 0; - if (nproc[2] > 0) inlet_layers_z = 0; + if (rank_info.ix > 0) inlet_layers_x = 0; + if (rank_info.jy > 0) inlet_layers_y = 0; + if (rank_info.kz > 0) inlet_layers_z = 0; + // outlet layers only apply to top part of domain + if (rank_info.ix < nproc[0]-1 ) outlet_layers_x = 0; + if (rank_info.jy < nproc[1]-1) outlet_layers_y = 0; + if (rank_info.kz < nproc[2]-1) outlet_layers_z = 0; // Fill remaining variables N = Nx*Ny*Nz; Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0; @@ -514,43 +525,46 @@ void Domain::CommInit() void Domain::ReadIDs(){ // Read the IDs from input file - int nprocs=nprocx()*nprocy()*nprocz(); - size_t readID; - char LocalRankString[8]; - char LocalRankFilename[40]; - //....................................................................... - if (rank() == 0) printf("Read input media... \n"); - //....................................................................... - sprintf(LocalRankString,"%05d",rank()); - sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); - // .......... READ THE INPUT FILE ....................................... - if (rank()==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n"); - sprintf(LocalRankFilename,"ID.%05i",rank()); - FILE *IDFILE = fopen(LocalRankFilename,"rb"); - if (IDFILE==NULL) ERROR("Domain::ReadIDs -- Error opening file: ID.xxxxx"); - readID=fread(id,1,N,IDFILE); - if (readID != size_t(N)) printf("Domain::ReadIDs -- Error reading ID (rank=%i) \n",rank()); - fclose(IDFILE); - // Compute the porosity - double sum; - double porosity; - double sum_local=0.0; - double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs); - if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6)); - //......................................................... - // If external boundary conditions are applied remove solid - if (BoundaryCondition > 0 && kproc() == 0){ - for (int k=0; k<3; k++){ - for (int j=0;j 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6)); + //......................................................... + // If external boundary conditions are applied remove solid + if (BoundaryCondition > 0 && kproc() == 0){ + if (inlet_layers_z < 4) inlet_layers_z=4; + for (int k=0; k 0 && kproc() == nprocz()-1){ - for (int k=Nz-3; k