Add naive Gray Scott simulation test

CMakeLists.txt

@@ -5,13 +5,21 @@ add_subdirectory(cmake)
 
 find_package(TBB COMPONENTS tbb REQUIRED)
 
+find_package(HDF5 COMPONENTS C CXX REQUIRED)
+
+include_directories(${HDF5_INCLUDE_DIRS})
+
+message(STATUS "HDF5_CXX_LIBRARIES = ${HDF5_CXX_LIBRARIES}")
+
 phoenix_base_project("PhoenixPerformance" "1.4.1"
 		"Set of performance test"
 		"https://gitlab.in2p3.fr/CTA-LAPP/PHOENIX_LIBS/PhoenixPerformance")
 
+pull_extra_module("OptionParser" "https://gitlab.in2p3.fr/CTA-LAPP/PHOENIX_LIBS/OptionParser.git")
 pull_extra_module("MicroBenchmark" "https://gitlab.in2p3.fr/CTA-LAPP/PHOENIX_LIBS/MicroBenchmark.git")
 pull_extra_module("TensorAlloc" "https://gitlab.in2p3.fr/CTA-LAPP/PHOENIX_LIBS/TensorAlloc.git")
 
+
 include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
 
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src ${TBB_TBB_INCLUDE_DIRS})

src/cpp20/Sequential/CMakeLists.txt

@@ -6,4 +6,6 @@ add_subdirectory(saxpy)
 add_subdirectory(barycenter)
 add_subdirectory(barycenter2d)
 
+add_subdirectory(GrayScott)
+
 

src/cpp20/Sequential/GrayScott/CMakeLists.txt

+project(Phoenix)
+cmake_minimum_required(VERSION 2.8)
+
+file(GLOB naiveSource "${CMAKE_CURRENT_SOURCE_DIR}/Naive/*.cpp")
+file(GLOB mainSource "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp")
+
+add_executable(naive_gray_scott ${naiveSource} ${mainSource})
+target_link_libraries(naive_gray_scott tensor_alloc option_parser data_stream string_utils ${HDF5_CXX_LIBRARIES} TBB::tbb)
+
+
+

src/cpp20/Sequential/GrayScott/MatrixHdf5.cpp

+/***************************************
+	Auteur : Pierre Aubert
+	Mail : aubertp7@gmail.com
+	Licence : CeCILL-C
+****************************************/
+
+
+//Warning : this file has been generated automatically by the phoenix_hdf5 program
+//You can find it at https://gitlab.in2p3.fr/CTA-LAPP/PHOENIX_LIBS/PhoenixHDF5
+//Do NOT modify it
+
+
+#include <string.h>
+#include <sstream>
+#include "MatrixHdf5.h"
+
+///Constructor of the class MatrixHdf5
+MatrixHdf5::MatrixHdf5(){
+	p__nbRow = 0lu;
+	p__tableName = "matrix";
+	p_nbCol = 0lu;
+	p_nbRow = 0lu;
+	p_tabMat = NULL;
+}
+
+///Destructor of the class MatrixHdf5
+MatrixHdf5::~MatrixHdf5(){
+	clear();
+}
+
+///Set the HDF5 name of the Table MatrixHdf5
+/**	@param name : name of the table to be saved
+*/void MatrixHdf5::setTableName(const std::string & name){
+	p__tableName = name;
+}
+
+///Get the total number of rows in the current Table MatrixHdf5
+/**	@return total number of rows
+*/size_t MatrixHdf5::getNbEntries() const{
+	return p__nbRow;
+}
+
+///Resize the table MatrixHdf5
+/*	@param nbRow : new number of rows of the Table
+*/
+void MatrixHdf5::resize(size_t nbRow){
+	if(nbRow == p__nbRow){return;}	//Nothing to do
+	clear();
+	allocate(nbRow);
+}
+
+///Clear the table MatrixHdf5 (delete all column)
+void MatrixHdf5::clear(){
+	if(p_tabMat != NULL){delete [] p_tabMat;p_tabMat = NULL;}
+}
+
+///Read the table MatrixHdf5 from given file
+/**	@param fileName : name of the HDF5 file to be used
+*/
+void MatrixHdf5::read(const std::string & fileName){
+	H5::H5File file(fileName, H5F_ACC_RDONLY);
+	read(file);
+}
+
+///Read the table MatrixHdf5 from given file
+/**	@param file : HDF5 file to be used
+*/
+void MatrixHdf5::read(const H5::H5File & file){
+	H5::DataSet dataset = file.openDataSet(p__tableName);
+	readDataSet(dataset);
+}
+
+///Read the table MatrixHdf5 from given group
+/**	@param group : HDF5 group to be used
+*/
+void MatrixHdf5::read(const H5::Group & group){
+	H5::DataSet dataset = group.openDataSet(p__tableName);
+	readDataSet(dataset);
+}
+
+///Create and write the table MatrixHdf5 in given file
+/**	@param fileName : name of the HDF5 file to be used
+*/
+void MatrixHdf5::write(const std::string & fileName) const{
+	H5::H5File file(fileName, H5F_ACC_TRUNC);
+	write(file);
+}
+
+///Create and write the table MatrixHdf5 in given file
+/**	@param file : HDF5 file to be used
+*/
+void MatrixHdf5::write(H5::H5File & file) const{
+	hsize_t dim[1];
+	dim[0] = p__nbRow;
+	H5::DataSpace space(1, dim);
+	H5::DataSet dataset = file.createDataSet(p__tableName, getCompTypeAll(), space);
+	writeDataSet(dataset);
+}
+
+///Create and write the table MatrixHdf5 in given file
+/**	@param group : HDF5 group to be used
+*/
+void MatrixHdf5::write(H5::Group & group) const{
+	hsize_t dim[1];
+	dim[0] = p__nbRow;
+	H5::DataSpace space(1, dim);
+	H5::DataSet dataset = group.createDataSet(p__tableName, getCompTypeAll(), space);
+	writeDataSet(dataset);
+}
+
+///Set a full row of the table MatrixHdf5
+/**	@param i : index of the row to be set
+ * 	@param tabMat : attribute to be set
+*/
+void MatrixHdf5::setRow(size_t i, float * tabMat){
+	setTabMat(i, tabMat);
+}
+
+///Get a full row of the table MatrixHdf5 (without tensor copy, only with pointer)
+/**	@param i : index of the row to get its values
+ * 	@param[out] tabMat : attribute to be get
+*/
+void MatrixHdf5::getRow(size_t i, float *& tabMat){
+	tabMat = getTabMat(i);
+}
+
+///Get a full row of the table MatrixHdf5 (without tensor copy, only with pointer)
+/**	@param i : index of the row to get its values
+ * 	@param[out] tabMat : attribute to be get
+*/
+void MatrixHdf5::getRow(size_t i, const float *& tabMat) const{
+	tabMat = getTabMat(i);
+}
+
+///Set the attribute tabMat (column image)
+/**	@param i : index of the row to be used
+ * 	@param tabVal : table of value to be copied
+*/
+void MatrixHdf5::setTabMat(size_t i, const float * tabVal){
+	size_t sizeRow(p_nbRow*p_nbCol);
+	memcpy(p_tabMat + i*sizeRow, tabVal, sizeRow);
+}
+
+///Get the full column of the attribute tabMat (column image)
+/**	@return pointer of the full column
+*/
+const float * MatrixHdf5::getTabMatFull() const{
+	return p_tabMat;
+}
+
+///Get the full column of the attribute tabMat (column image)
+/**	@return pointer of the full column
+*/
+float * MatrixHdf5::getTabMatFull(){
+	return p_tabMat;
+}
+
+///Get the tensor i of the attribute tabMat (column image)
+/**	@param i : index of the row to be used
+ * 	@return pointer to the corresponding tensor
+*/
+const float * MatrixHdf5::getTabMat(size_t i) const{
+	return p_tabMat + i*p_nbRow*p_nbCol;
+}
+
+///Get the tensor i of the attribute tabMat (column image)
+/**	@param i : index of the row to be used
+ * 	@return pointer to the corresponding tensor
+*/
+float * MatrixHdf5::getTabMat(size_t i){
+	return p_tabMat + i*p_nbRow*p_nbCol;
+}
+
+///Set all the dimentions of the Tensor in the table
+/**	@param nbCol : set the tensor dimention nbCol
+ * 	@param nbRow : set the tensor dimention nbRow
+*/
+void MatrixHdf5::setAllDim(size_t nbCol, size_t nbRow){
+	setNbCol(nbCol);
+	setNbRow(nbRow);
+}
+
+///Set the Tensor dimention nbCol
+/**	@param val : set the tensor dimention nbCol
+*/
+void MatrixHdf5::setNbCol(size_t val){
+	 p_nbCol = val;
+}
+
+///Get the Tensor dimention nbCol
+/**	@return the tensor dimention nbCol
+*/
+size_t MatrixHdf5::getNbCol() const{
+	return p_nbCol;
+}
+
+///Set the Tensor dimention nbRow
+/**	@param val : set the tensor dimention nbRow
+*/
+void MatrixHdf5::setNbRow(size_t val){
+	 p_nbRow = val;
+}
+
+///Get the Tensor dimention nbRow
+/**	@return the tensor dimention nbRow
+*/
+size_t MatrixHdf5::getNbRow() const{
+	return p_nbRow;
+}
+
+///Get the offset of the attribute tabMat
+/**	@return offset of the attribute in bytes
+*/
+size_t MatrixHdf5::getOffsetTabMat() const{
+	return 0lu;
+}
+
+H5::CompType MatrixHdf5::getCompTypeAll() const{
+	size_t sizeAll(sizeof(float)*p_nbRow*p_nbCol);
+	H5::CompType typeCol(sizeAll);
+	typeCol.insertMember("image", getOffsetTabMat(), getTypeTabMat());
+	return typeCol;
+}
+
+///Get DataType of attribute tabMat (column image)
+/**	@return DataType of the attribute tabMat
+*/
+H5::CompType MatrixHdf5::getCompTypeTabMat() const{
+	H5::CompType typeCol(sizeof(float)*p_nbRow*p_nbCol);
+	typeCol.insertMember("image", 0, getTypeTabMat());
+	return typeCol;
+}
+
+///Get DataType of attribute tabMat (column image)
+/**	@return DataType of the attribute tabMat
+*/
+H5::DataType MatrixHdf5::getTypeTabMat() const{
+	hsize_t dims[2];
+	dims[0] = p_nbRow;
+	dims[1] = p_nbCol;
+	H5::ArrayType arrayType(H5::PredType::NATIVE_FLOAT, 2, dims);
+	return arrayType;
+}
+
+///Read the given DataSet and fill the Table with it
+/**	@param dataset : dataset to be used
+*/
+void MatrixHdf5::readDataSet(const H5::DataSet & dataset){
+	H5::CompType compType = dataset.getCompType();
+	readDimTabMat(compType);
+	H5::DataSpace dataSpace = dataset.getSpace();
+	size_t nbEntries(dataSpace.getSimpleExtentNpoints());
+	resize(nbEntries);
+	dataset.read(p_tabMat, getCompTypeTabMat());
+}
+
+///Write the given DataSet and fill the Table with it
+/**	@param[out] dataset : dataset to be modified
+*/
+void MatrixHdf5::writeDataSet(H5::DataSet & dataset) const{
+	dataset.write(p_tabMat, getCompTypeTabMat());
+}
+
+///Allocate the table MatrixHdf5 (delete all column)
+/*	@param nbRow : new number of rows of the Table
+*/
+void MatrixHdf5::allocate(size_t nbRow){
+	p__nbRow = nbRow;
+	p_tabMat = new float[p__nbRow*p_nbRow*p_nbCol];
+}
+
+///Update the dimentions of the tensor tabMat (column 'image')
+/**	@param compType : main type to be used to update the tensor dimentions*/
+void MatrixHdf5::readDimTabMat(const H5::CompType & compType){
+	int indexCol = compType.getMemberIndex("image");
+	H5::ArrayType arrayType = compType.getMemberArrayType(indexCol);
+	//Check if the number of dientions matches
+	if(arrayType.getArrayNDims() != 2){
+		std::stringstream strError;
+		strError << "MatrixHdf5::readDimTabMat wrong number of dimentions : expect 2 but "<< arrayType.getArrayNDims() <<" provided from the file" << std::endl;
+		throw std::runtime_error(strError.str());
+	}
+	hsize_t dims[2];
+	arrayType.getArrayDims(dims);
+	p_nbRow = dims[0];
+	p_nbCol = dims[1];
+}
+

src/cpp20/Sequential/GrayScott/MatrixHdf5.h

+/***************************************
+	Auteur : Pierre Aubert
+	Mail : aubertp7@gmail.com
+	Licence : CeCILL-C
+****************************************/
+
+
+//Warning : this file has been generated automatically by the phoenix_hdf5 program
+//You can find it at https://gitlab.in2p3.fr/CTA-LAPP/PHOENIX_LIBS/PhoenixHDF5
+//Do NOT modify it
+
+
+#ifndef __MATRIXHDF5_H__
+#define __MATRIXHDF5_H__
+
+#include <iostream>
+#include "H5Cpp.h"
+
+///@brief Matrix of the concentation of the U specie
+class MatrixHdf5{
+	public:
+		MatrixHdf5();
+		virtual ~MatrixHdf5();
+
+		void setTableName(const std::string & name);
+		size_t getNbEntries() const;
+		void resize(size_t nbRow);
+		void clear();
+
+		void read(const std::string & fileName);
+		void read(const H5::H5File & file);
+		void read(const H5::Group & group);
+
+		void write(const std::string & fileName) const;
+		void write(H5::H5File & file) const;
+		void write(H5::Group & group) const;
+
+		void setRow(size_t i, float * tabMat);
+		void getRow(size_t i, float *& tabMat);
+		void getRow(size_t i, const float *& tabMat) const;
+		void setTabMat(size_t i, const float * tabVal);
+		const float * getTabMatFull() const;
+		float * getTabMatFull();
+		const float * getTabMat(size_t i) const;
+		float * getTabMat(size_t i);
+
+		size_t getOffsetTabMat() const;
+		H5::CompType getCompTypeAll() const;
+		H5::CompType getCompTypeTabMat() const;
+		H5::DataType getTypeTabMat() const;
+
+		void setAllDim(size_t nbCol, size_t nbRow);
+		void setNbCol(size_t val);
+		size_t getNbCol() const;
+		void setNbRow(size_t val);
+		size_t getNbRow() const;
+
+	private:
+		void readDataSet(const H5::DataSet & dataset);
+		void writeDataSet(H5::DataSet & dataset) const;
+		void allocate(size_t nbRow);
+
+		void readDimTabMat(const H5::CompType & compType);
+
+		///Number of rows in the table MatrixHdf5
+		size_t p__nbRow;
+
+		///HDF5 name of the table MatrixHdf5
+		std::string p__tableName;
+
+		///Tensor dimension nbCol
+		size_t p_nbCol;
+		///Tensor dimension nbRow
+		size_t p_nbRow;
+
+		///Matrix of the concentration of the U specie
+		float * p_tabMat;
+};
+
+#endif
+

src/cpp20/Sequential/GrayScott/Naive/naive_propagation.cpp

@@ -4,7 +4,7 @@
 	Licence : CeCILL-C
 ****************************************/
 
-
+#include <algorithm>
 #include "naive_propagation.h"
 
 
@@ -21,14 +21,39 @@
  * 	@param diffudionRateU : diffusion rate of the U specie
  * 	@param diffudionRateV : diffusion rate of the V specie
  * 	@param feedRate : rate of the process which feeds U and drains U, V and P
- * 	@param killRate : rate of hte process which converts V into P
+ * 	@param killRate : rate of the process which converts V into P
 */
-void naive_propagation(float * outMatU, float * outMatV, const float * matU, const float * matV, size_t nbRow, size_t nbCol,
-		       const float * matDeltaSquare, size_t nbStencilRow, size_t nbStencilCol,
+void naive_propagation(float * outMatU, float * outMatV, const float * matU, const float * matV, long nbRow, long nbCol,
+		       const float * matDeltaSquare, long nbStencilRow, long nbStencilCol,
 		       float diffudionRateU, float diffusionRateV, float feedRate, float killRate)
 {
+	long offsetStencilRow((nbStencilRow - 1l)/2l);
+	long offsetStencilCol((nbStencilCol - 1l)/2l);
 	
-	
-	
+	for(long i(0l); i < nbRow; ++i){
+		for(long j(0l); j < nbCol; ++j){
+			float fullU(0.0f), fullV(0.0f);
+			
+			long firstRowStencil(std::max(i - offsetStencilRow, 0l));
+			long firstColStencil(std::max(j - offsetStencilCol, 0l));
+			
+			long lastRowStencil(std::min(i + offsetStencilRow + 1l, nbStencilRow));
+			long lastColStencil(std::min(j + offsetStencilCol + 1l, nbStencilCol));
+			
+			for(long k(firstRowStencil); k < lastRowStencil; ++k){
+				for(long l(firstColStencil); l < lastColStencil; ++l){
+					fullU += matU[k*nbCol + l];
+					fullV += matV[k*nbCol + l];
+				}
+			}
+			float u(matU[i*nbCol + j]), v(matV[i*nbCol + j]);
+			float uvSquare(u*v*v);
+			float du(diffudionRateU*fullU - uvSquare + feedRate*(1 - u));
+			float dv(diffusionRateV*fullV + uvSquare - (feedRate + killRate)*v);
+			
+			outMatU[i*nbCol + j] = u + du;
+			outMatV[i*nbCol + j] = v + dv;
+		}
+	}
 }
 

src/cpp20/Sequential/GrayScott/Naive/naive_propagation.h

@@ -10,8 +10,8 @@
 #include <iostream>
 
 
-void naive_propagation(float * outMatU, float * outMatV, const float * matU, const float * matV, size_t nbRow, size_t nbCol,
-		       const float * matDeltaSquare, size_t nbStencilRow, size_t nbStencilCol,
+void naive_propagation(float * outMatU, float * outMatV, const float * matU, const float * matV, long nbRow, long nbCol,
+		       const float * matDeltaSquare, long nbStencilRow, long nbStencilCol,
 		       float diffudionRateU, float diffusionRateV, float feedRate, float killRate);
 
 #endif

src/cpp20/Sequential/GrayScott/main.cpp

+/***************************************
+	Auteur : Pierre Aubert
+	Mail : aubertp7@gmail.com
+	Licence : CeCILL-C
+****************************************/
+
+#include "OptionParser.h"
+#include "PTensor.h"
+
+#include "Naive/naive_propagation.h"
+#include "MatrixHdf5.h"
+
+///Create the OptionParser of this program
+/**	@return OptionParser of this program
+*/
+OptionParser createOptionParser(){
+	OptionParser parser(true, __PROGRAM_VERSION__);
+	parser.setExampleLongOption("naive_gray_scott --killrate=0.062 --feedrate=0.03 --nbimage=100 --nbrow=1080 --nbcol=1920 --output=outputFile.hdf5");
+	parser.setExampleShortOption("naive_gray_scott -k 0.062 -f 0.03 -n 100 -r 1080 -c 1920 -o outputFile.hdf5");
+	
+	float killRate(0.062f), feedRate(0.03f);
+	size_t nbImage(100lu), nbRow(1080lu), nbCol(1920lu);
+	parser.addOption("killrate", "k", killRate, "rate of the process which converts V into P");
+	parser.addOption("feedrate", "f", feedRate, "rate of the process which feeds U and drains U, V and P");
+	parser.addOption("nbimage", "n", nbImage, "number of images to be created");
+	
+	parser.addOption("nbrow", "r", nbRow, "number of rows of the images to be created");
+	parser.addOption("nbcol", "c", nbCol, "number of columns of the images to be created");
+	
+	std::string defaultOutputFile("output.h5");
+	parser.addOption("output", "o", defaultOutputFile, "Output file to be created with results");
+	
+	return parser;
+}
+
+///Swap two values
+/**	@param[out] a : value will be b
+ * 	@param[out] b : value will be a
+*/
+template<typename T>
+void swapValue(T & a, T & b){
+	T c(a);
+	a = b;
+	b = c;
+}
+
+///Simulate the images
+/**	@param nbRow : number of rows of the images to be created
+ * 	@param nbCol : number of columns of the images to be created
+ * 	@param nbImage : number of images to be created
+ * 	@param killRate : rate of the process which converts V into P
+ * 	@param feedRate : rate of the process which feeds U and drains U, V and P
+ * 	@param outputFile : name of the file to be created
+ * 	@return true on succsess, false otherwise
+*/
+bool simulateImage(size_t nbRow, size_t nbCol, size_t nbImage, float killRate, float feedRate, const std::string & outputFile){
+	
+	MatrixHdf5 fullMat;
+	fullMat.setAllDim(nbCol, nbRow);
+	fullMat.resize(nbImage);
+	
+	PTensor<float> tmpInU(AllocMode::ALIGNED, nbRow, nbCol);
+	PTensor<float> tmpInV(AllocMode::ALIGNED, nbRow, nbCol);
+	PTensor<float> tmpOutU(AllocMode::ALIGNED, nbRow, nbCol);
+	PTensor<float> tmpOutV(AllocMode::ALIGNED, nbRow, nbCol);
+	tmpInU.fill(0.0f);
+	tmpOutU.fill(0.0f);
+	tmpInV.fill(1.0f);
+	tmpOutV.fill(1.0f);
+	
+	float *tmpU1 = tmpInU.getData(), *tmpU2 = tmpOutU.getData();
+	float *tmpV1 = tmpInV.getData(), *tmpV2 = tmpOutV.getData();
+	
+	long nbStencilRow(3l), nbStencilCol(3l);
+	
+	float diffudionRateU(0.05f), diffusionRateV(0.05f);
+	float matDeltaSquare[] = 	{0.05f, 0.2f, 0.05f,
+					0.2f, -1.0f, 0.2f,
+					0.05f, 0.2f, 0.05f};
+	
+	for(size_t i(0lu); i < nbImage; ++i){
+		std::cout << "simulateImage n°" << i << "..." << std::endl;
+		float * matInputU = tmpU2;
+		float * matInputV = tmpV2;
+		float * matOutputU = tmpU1;
+		float * matOutputV = tmpV1;
+		if(i%2lu == 0lu){
+			matInputU = tmpU1;
+			matInputV = tmpV1;
+			matOutputU = tmpU2;
+			matOutputV = tmpV2;
+		}
+		naive_propagation(matOutputU, matOutputV, matInputU, matInputV, nbRow, nbCol,
+				matDeltaSquare, nbStencilRow, nbStencilCol,
+				diffudionRateU, diffusionRateV, feedRate, killRate);
+		
+		fullMat.setRow(i, matOutputU);
+		
+		///Let's swap the pointer
+		swapValue(tmpU1, tmpU2);
+		swapValue(tmpV1, tmpV2);
+	}
+	std::cerr << "Done" << std::endl;
+	//Let's save the output file
+	fullMat.write(outputFile);
+	return true;
+}
+
+
+int main(int argc, char** argv){
+	OptionParser parser = createOptionParser();
+	parser.parseArgument(argc, argv);
+	
+	const OptionMode & defaultMode = parser.getDefaultMode();
+	float killRate(0.062f), feedRate(0.03f);
+	size_t nbImage(100lu), nbRow(1080lu), nbCol(1920lu);
+	defaultMode.getValue(killRate, "killrate");
+	defaultMode.getValue(feedRate, "feedrate");
+	defaultMode.getValue(nbImage, "nbimage");