Still some trouble with the Gray Scott simulation

src/cpp20/Sequential/GrayScott/MatrixHdf5.cpp

@@ -138,7 +138,7 @@ void MatrixHdf5::getRow(size_t i, const float *& tabMat) const{
 */
 void MatrixHdf5::setTabMat(size_t i, const float * tabVal){
 	size_t sizeRow(p_nbRow*p_nbCol);
-	memcpy(p_tabMat + i*sizeRow, tabVal, sizeRow);
+	memcpy(p_tabMat + i*sizeRow, tabVal, sizeRow*sizeof(float));
 }
 
 ///Get the full column of the attribute tabMat (column image)

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

@@ -22,10 +22,11 @@
  * 	@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 the process which converts V into P
+ * 	@param dt : time interval between two steps
 */
 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)
+		       float diffudionRateU, float diffusionRateV, float feedRate, float killRate, float dt)
 {
 	long offsetStencilRow((nbStencilRow - 1l)/2l);
 	long offsetStencilCol((nbStencilCol - 1l)/2l);
@@ -40,19 +41,24 @@ void naive_propagation(float * outMatU, float * outMatV, const float * matU, con
 			long lastRowStencil(std::min(i + offsetStencilRow + 1l, nbStencilRow));
 			long lastColStencil(std::min(j + offsetStencilCol + 1l, nbStencilCol));
 			
+			long stencilIndexRow(0l), stencilIndexCol(0l);
+			float u(matU[i*nbCol + j]), v(matV[i*nbCol + j]);
 			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 deltaSquare(matDeltaSquare[stencilIndexRow*nbStencilCol + stencilIndexCol]);
+					float deltaSquare(1.0f);
+					fullU += (matU[k*nbCol + l] - u)*deltaSquare;
+					fullV += (matV[k*nbCol + l] - v)*deltaSquare;
+					++stencilIndexCol;
 				}
+				++stencilIndexRow;
 			}
-			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);
+			float du(diffudionRateU*fullU*u - uvSquare + feedRate*(1 - u));
+			float dv(diffusionRateV*fullV*v + uvSquare - (feedRate + killRate)*v);
 			
-			outMatU[i*nbCol + j] = u + du;
-			outMatV[i*nbCol + j] = v + dv;
+			outMatU[i*nbCol + j] = u + du*dt;
+			outMatV[i*nbCol + j] = v + dv*dt;
 		}
 	}
 }

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

@@ -12,6 +12,6 @@
 
 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);
+		       float diffudionRateU, float diffusionRateV, float feedRate, float killRate, float dt);
 
 #endif

src/cpp20/Sequential/GrayScott/main.cpp

@@ -18,18 +18,25 @@ OptionParser createOptionParser(){
 	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);
+	float killRate(0.054f), feedRate(0.014f);
 	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");
+	size_t nbExtraStep(1lu);
+	parser.addOption("nbextrastep", "e", nbExtraStep, "number of extra steps to be computed between images");
 	
 	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");
+	float dt(1.0f);
+	parser.addOption("deltat", "t", dt, "time interval between two computation");
 	
 	std::string defaultOutputFile("output.h5");
 	parser.addOption("output", "o", defaultOutputFile, "Output file to be created with results");
 	
+// 	0.014,0.054,
+// 	0,#000000,0.2,#00FF00,0.21,#FFFF00,0.4,#FF0000,0.6,#FFFFFF
+	
 	return parser;
 }
 
@@ -48,12 +55,15 @@ void swapValue(T & a, T & b){
 /**	@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 nbExtraStep : number of extra steps to be computed between images
  * 	@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 dt : time interval between two computation
  * 	@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){
+bool simulateImage(size_t nbRow, size_t nbCol, size_t nbImage, size_t nbExtraStep, float killRate, float feedRate, float dt, const std::string & outputFile){
+	std::cout << "simulateImage : nbRow = " << nbRow << ", nbCol = " << nbCol << std::endl;
 	
 	MatrixHdf5 fullMat;
 	fullMat.setAllDim(nbCol, nbRow);
@@ -63,48 +73,57 @@ bool simulateImage(size_t nbRow, size_t nbCol, size_t nbImage, float killRate, f
 	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);
+	tmpInU.fill(1.0f);
+	tmpOutU.fill(1.0f);
+	tmpInV.fill(0.0f);
+	tmpOutV.fill(0.0f);
 	
-	tmpInU.setValue(nbRow/2lu, nbCol/2lu, 1.0f);
-	tmpInU.setValue(nbRow/3lu, nbCol/2lu, 1.0f);
-	tmpInU.setValue(nbRow/2lu, nbCol/3lu, 1.0f);
-	tmpInU.setValue(nbRow/4lu, nbCol/2lu, 1.0f);
-	tmpInU.setValue(nbRow/2lu, nbCol/4lu, 1.0f);
+	for(size_t i((2lu*nbRow)/5lu); i < (3lu*nbRow)/5lu; ++i){
+		for(size_t j((2lu*nbCol)/5lu); j < (3lu*nbCol)/5lu; ++j){
+			tmpInU.setValue(i, j, 0.0f);
+			tmpInV.setValue(i, j, 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};
-	
+	float diffudionRateU(0.01f), diffusionRateV(0.5f);
+// 	float matDeltaSquare[] = 	{0.05f, 0.2f, 0.05f,
+// 					0.2f, -1.0f, 0.2f,
+// 					0.05f, 0.2f, 0.05f};
+	float matDeltaSquare[] = 	{1.0f, 1.0f, 1.0f,
+					1.0f, 1.0f, 1.0f,
+					1.0f, 1.0f, 1.0f};
+					
 	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;
+		for(size_t j(0lu); j < nbExtraStep; ++j){
+			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, dt);
+			
+			fullMat.setRow(i, matOutputV);
+			
+			///Let's swap the pointer
+	// 		swapValue(tmpU1, tmpU2);
+	// 		swapValue(tmpV1, tmpV2);
+			
+			tmpInU = tmpOutU;
+			tmpInV = tmpOutV;
 		}
-		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
@@ -118,18 +137,21 @@ int main(int argc, char** argv){
 	parser.parseArgument(argc, argv);
 	
 	const OptionMode & defaultMode = parser.getDefaultMode();
-	float killRate(0.062f), feedRate(0.03f);
-	size_t nbImage(100lu), nbRow(1080lu), nbCol(1920lu);
+	float killRate(0.062f), feedRate(0.03f), dt(1.0f);
+	size_t nbImage(100lu), nbRow(1080lu), nbCol(1920lu), nbExtraStep(1lu);
 	defaultMode.getValue(killRate, "killrate");
 	defaultMode.getValue(feedRate, "feedrate");
 	defaultMode.getValue(nbImage, "nbimage");
+	defaultMode.getValue(nbExtraStep, "nbextrastep");
 	defaultMode.getValue(nbRow, "nbrow");
 	defaultMode.getValue(nbCol, "nbcol");
+	defaultMode.getValue(dt, "deltat");
+	
 	
 	std::string outputFile("");
 	defaultMode.getValue(outputFile, "output");
 	
-	bool b(simulateImage(nbRow, nbCol, nbImage, killRate, feedRate, outputFile));
+	bool b(simulateImage(nbRow, nbCol, nbImage, nbExtraStep, killRate, feedRate, dt, outputFile));
 	return b - 1;
 }