lbpm_nonnewtonian_simulator.h

// 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 1000






enum AnalysisType{ AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02, 
	CopySimState=0x04, ComputeAverages=0x08, CreateRestart=0x10, WriteVis=0x20 };








template<class TYPE>
void DeleteArray( const TYPE *p )
{
	delete [] p;
}








// 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
{
public:
	WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cDen_,
			std::shared_ptr<double> cDistEven_, std::shared_ptr<double>cDistOdd_, int N_ ):
				filename(filename_), cDen(cDen_), cDistEven(cDistEven_), cDistOdd(cDistOdd_), N(N_) {}
	virtual void run() {
		PROFILE_START("Save Checkpoint",1);
		WriteCheckpoint(filename,cDen.get(),cDistEven.get(),cDistOdd.get(),N);
		PROFILE_STOP("Save Checkpoint",1);
	};
	virtual bool has_result() const { return false; }
private:
	WriteRestartWorkItem();
	const char* filename;
	std::shared_ptr<double> cDen, cDistEven, cDistOdd;
	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;
typedef std::shared_ptr<std::vector<BlobIDType> > BlobIDList;



//
//class BlobIdentificationWorkItem1: public ThreadPool::WorkItem
//{
//public:
//    BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
//        std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
//        BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ):
//        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_)
//        {
//            MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
//        }
//    ~BlobIdentificationWorkItem1() { MPI_Comm_free(&newcomm); }
//    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,newcomm);
//        PROFILE_STOP("Identify blobs",1);
//    }
//    virtual bool has_result() const { return false; }
//private:
//    BlobIdentificationWorkItem1();
//    int timestep;
//    int Nx, Ny, Nz;
//    const RankInfoStruct& rank_info;
//    std::shared_ptr<const DoubleArray> phase;
//    const DoubleArray& dist;
//    BlobIDstruct last_id, new_index, new_id;
//    BlobIDList new_list;
//    MPI_Comm newcomm;
//};
//


//
//
//class BlobIdentificationWorkItem2: public ThreadPool::WorkItem
//{
//public:
//    BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_,
//        std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
//        BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ ):
//        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_)
//        {
//            MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
//        }
//    ~BlobIdentificationWorkItem2() { MPI_Comm_free(&newcomm); }
//    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,newcomm);
//            // 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);
//    }
//    virtual bool has_result() const { return false; }
//private:
//    BlobIdentificationWorkItem2();
//    int timestep;
//    int Nx, Ny, Nz;
//    const RankInfoStruct& rank_info;
//    std::shared_ptr<const DoubleArray> phase;
//    const DoubleArray& dist;
//    BlobIDstruct last_id, new_index, new_id;
//    BlobIDList new_list;
//    MPI_Comm newcomm;
//};
//


// Helper class to write the vis file from a thread
class WriteVisWorkItem: public ThreadPool::WorkItem
{
public:
	WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
			TwoPhase& Avgerages_, fillHalo<double>& fillData_ ):
				timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_)
{
		MPI_Comm_dup(MPI_COMM_WORLD,&newcomm);
}
	~WriteVisWorkItem() { MPI_Comm_free(&newcomm); }
	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<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);
		IO::writeData( timestep, visData, newcomm );
		PROFILE_STOP("Save Vis",1);
	};
	virtual bool has_result() const { return false; }
private:
	WriteVisWorkItem();
	int timestep;
	std::vector<IO::MeshDataStruct>& visData;
	TwoPhase& Averages;
	fillHalo<double>& fillData;
	MPI_Comm newcomm;
};



// Helper class to run the analysis from within a thread
// Note: Averages will be modified after the constructor is called
class AnalysisWorkItem: public ThreadPool::WorkItem
{
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 ( (type&CopyPhaseIndicator) != 0 ) {
			// Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
		}
		if ( (type&ComputeAverages) != 0 ) {
			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);
		}
	}
	virtual bool has_result() const { return false; }
private:
	AnalysisWorkItem();
	AnalysisType type;
	int timestep;
	TwoPhase& Averages;
	BlobIDstruct blob_ids;
	BlobIDList id_list;
	double beta;
};


// Function to start the analysis
void run_analysis( int timestep, int restart_interval, 
		const RankInfoStruct& rank_info, TwoPhase& Averages,
		BlobIDstruct& last_ids, BlobIDstruct& last_index, BlobIDList& last_id_map,
		int Nx, int Ny, int Nz, bool pBC, double err,
		double *Pressure, const double *Velocity,
		const char *ID, const double *f_even, const double *f_odd,
		const char *LocalRestartFile, std::vector<IO::MeshDataStruct>& visData, fillHalo<double>& fillData,
		ThreadPool& tpool, AnalysisWaitIdStruct& wait )
{

	int N = Nx*Ny*Nz;

	/*
	 * 2
	 */

	// Determin the analysis we want to perform
	AnalysisType type = AnalyzeNone;
	//    if ( timestep%ANALYSIS_INTERVAL + 5 == ANALYSIS_INTERVAL ) {
	//        // Copy the phase indicator field for the earlier timestep
	//        type = static_cast<AnalysisType>( type | CopyPhaseIndicator );
	//    }
	//    if ( timestep%BLOBID_INTERVAL == 0 ) {
	//        // Identify blobs and update global ids in time
	//        type = static_cast<AnalysisType>( type | IdentifyBlobs );
	//    }
	//
	//    if ( timestep%ANALYSIS_INTERVAL == 0 ) {
	//        // Copy the averages to the CPU (and identify blobs)
	//        type = static_cast<AnalysisType>( type | CopySimState );
	//        type = static_cast<AnalysisType>( type | IdentifyBlobs );
	//    }
	//    if ( timestep%ANALYSIS_INTERVAL == 5 ) {
	//        // Run the analysis
	//        type = static_cast<AnalysisType>( type | ComputeAverages );
	//    }
	//    if (timestep%restart_interval == 0) {
	//        // 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 | CopySimState );
		type = static_cast<AnalysisType>( type | IdentifyBlobs );
	}

	// Return if we are not doing anything
	if ( type == AnalyzeNone )
		return;

	/*
	 * 3
	 */

	PROFILE_START("start_analysis");

	// Copy the appropriate variables to the host (so we can spawn new threads)
	ScaLBL_DeviceBarrier();

	/*
	 * 4
	 */

	PROFILE_START("Copy data to host",1);

	std::shared_ptr<DoubleArray> phase;
	//    if ( (type&CopyPhaseIndicator)!=0 || (type&ComputeAverages)!=0 ||  (type&CopySimState)!=0 || (type&IdentifyBlobs)!=0 )
	//    {
	//        phase = std::shared_ptr<DoubleArray>(new DoubleArray(Nx,Ny,Nz));
	//        ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double));
	//    }


	//    if ( (type&CopyPhaseIndicator)!=0 ) {
	//        memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double));
	//        //Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
	//    }


	//    if ( (type&ComputeAverages)!=0 ) {
	//        memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double));
	//        //Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tminus);
	//    }


	if ( (type&CopySimState) != 0 ) {
		// Copy the members of Averages to the cpu (phase was copied above)
		// Wait
		PROFILE_START("Copy-Pressure",1);
		ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
		ScaLBL_DeviceBarrier();
		PROFILE_STOP("Copy-Pressure",1);
		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-State",1);
		//  memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
		ScaLBL_CopyToHost(Averages.Press.data(),Pressure,N*sizeof(double));
		ScaLBL_CopyToHost(Averages.Vel_x.data(),&Velocity[0],N*sizeof(double));
		ScaLBL_CopyToHost(Averages.Vel_y.data(),&Velocity[N],N*sizeof(double));
		ScaLBL_CopyToHost(Averages.Vel_z.data(),&Velocity[2*N],N*sizeof(double));
		PROFILE_STOP("Copy-State",1);
	}

	std::shared_ptr<double> cDen, cDistEven, cDistOdd;
	if ( (type&CreateRestart) != 0 ) {
		// Copy restart data to the CPU
		cDen = std::shared_ptr<double>(new double[2*N],DeleteArray<double>);
		cDistEven = std::shared_ptr<double>(new double[10*N],DeleteArray<double>);
		cDistOdd = std::shared_ptr<double>(new double[9*N],DeleteArray<double>);
		ScaLBL_CopyToHost(cDistEven.get(),f_even,10*N*sizeof(double));
		ScaLBL_CopyToHost(cDistOdd.get(),f_odd,9*N*sizeof(double));
		//        ScaLBL_CopyToHost(cDen.get(),Den,2*N*sizeof(double));
	}


	PROFILE_STOP("Copy data to host",1);


	/*
	 * 5
	 */

	//
	//    // Spawn threads to do blob identification work
	//    if ( (type&IdentifyBlobs)!=0 ) {
	//        BlobIDstruct new_index(new std::pair<int,IntArray>(0,IntArray()));
	//        BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray()));
	//        BlobIDList new_list(new std::vector<BlobIDType>());
	//        ThreadPool::WorkItem *work1 = new BlobIdentificationWorkItem1(timestep,
	//            Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_index,new_ids,new_list);
	//        ThreadPool::WorkItem *work2 = new BlobIdentificationWorkItem2(timestep,
	//            Nx,Ny,Nz,rank_info,phase,Averages.SDs,last_ids,new_index,new_ids,new_list);
	//        work1->add_dependency(wait.blobID);
	//        work2->add_dependency(tpool.add_work(work1));
	//        wait.blobID = tpool.add_work(work2);
	//        last_index = new_index;
	//        last_ids = new_ids;
	//        last_id_map = new_list;
	//    }
	//

	/*
	 * 6
	 */

	//    // Spawn threads to do the analysis work
	//    if ( (type&ComputeAverages) != 0 ) {
	//    	double beta = 0.0;
	//        ThreadPool::WorkItem *work = new AnalysisWorkItem(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);
	//    }
	//

	/*
	 * 7
	 */

	// Spawn a thread to write the restart file
	if ( (type&CreateRestart) != 0 ) {
		int rank = MPI_WORLD_RANK();

		// Wait for previous restart files to finish writing (not necessary, but helps to ensure memory usage is limited)
		tpool.wait(wait.restart);
		// Retain the timestep associated with the restart files
		if (rank==0){
			FILE *Rst = fopen("Restart.txt","w");
			fprintf(Rst,"%i\n",timestep+5);
			fclose(Rst);
		}
		// 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);
	}


	/*
	 * 8
	 */

	// 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 );
		WriteVisWorkItem *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");
	printf("Returning from the run_analysis function...\n");
}




//if (pBC) {
//err = fabs(sat_w - sat_w_previous);
//sat_w_previous = sat_w;
//if (rank==0){
// printf("Timestep %i: change in saturation since last checkpoint is %f \n",timestep,err);
// }
// }








/*    #ifdef USE_CUDA
    if ( tpool.getQueueSize()<=3 && tpool.getNumThreads()>0 && 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
        // Add more detailed estimates of the update frequency required to track blobs
        type = static_cast<AnalysisType>( type | IdentifyBlobs );
    }
#endif
 */