Commit 82250241 authored by Pierre Aubert's avatar Pierre Aubert
Browse files

Seems to be bugs in the intrinsics version of the Gray Scott reaction

parent ee30d2f9
...@@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 2.8) ...@@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 2.8)
set(CONFIG_GRAYSCOTT "5, 7, 10, 20, 30, 40, 50, 60") set(CONFIG_GRAYSCOTT "5, 7, 10, 20, 30, 40, 50, 60")
set(EXTRA_DEPENDENCIES gray_scott_data_format tensor_alloc data_stream TBB::tbb) set(EXTRA_DEPENDENCIES gray_scott_data_format tensor_alloc data_stream ${HDF5_CXX_LIBRARIES} TBB::tbb)
set(progNaiveSrc ${CMAKE_CURRENT_SOURCE_DIR}/../Naive/naive_propagation.cpp main.cpp) set(progNaiveSrc ${CMAKE_CURRENT_SOURCE_DIR}/../Naive/naive_propagation.cpp main.cpp)
......
...@@ -46,8 +46,10 @@ void evaluateGrayScott(size_t nbElement){ ...@@ -46,8 +46,10 @@ void evaluateGrayScott(size_t nbElement){
tmpV2 = tmpVecOutV.getData(); tmpV2 = tmpVecOutV.getData();
float * ptrVecMatStencil = vecMatDeltaSquare.getData(); float * ptrVecMatStencil = vecMatDeltaSquare.getData();
size_t nbVecRow(tmpVecInV.getFullNbRow()), nbVecCol(tmpVecInV.getNbCol());
micro_benchmarkAutoNsPrint("evaluate GrayScott reaction, intrinsics", nbElement, grayscott_propagation, micro_benchmarkAutoNsPrint("evaluate GrayScott reaction, intrinsics", nbElement, grayscott_propagation,
tmpU2, tmpV2, tmpU1, tmpV1, nbRow, nbCol, 0l, tmpU2, tmpV2, tmpU1, tmpV1, nbVecRow, nbVecCol, 0l,
ptrVecMatStencil, nbStencilRow, nbStencilCol, ptrVecMatStencil, nbStencilRow, nbStencilCol,
diffudionRateU, diffusionRateV, feedRate, killRate, dt); diffudionRateU, diffusionRateV, feedRate, killRate, dt);
} }
......
...@@ -98,10 +98,10 @@ void grayscott_propagation(float * outMatVecU, float * outMatVecV, const float * ...@@ -98,10 +98,10 @@ void grayscott_propagation(float * outMatVecU, float * outMatVecV, const float *
PRegVecf vecDvDt(plib_mul_ps(vecDv, vecDt)); PRegVecf vecDvDt(plib_mul_ps(vecDv, vecDt));
PRegVecf vecUPlusDuDt(plib_add_ps(vecU, vecDuDt)); PRegVecf vecUPlusDuDt(plib_add_ps(vecU, vecDuDt));
PRegVecf vecVPluDvDt(plib_add_ps(vecV, vecDvDt)); PRegVecf vecVPlusDvDt(plib_add_ps(vecV, vecDvDt));
plib_store_ps(outMatVecU + (i*nbVecCol + j)*PLIB_VECTOR_SIZE_FLOAT, vecUPlusDuDt); plib_store_ps(outMatVecU + (i*nbVecCol + j)*PLIB_VECTOR_SIZE_FLOAT, vecUPlusDuDt);
plib_store_ps(outMatVecV + (i*nbVecCol + j)*PLIB_VECTOR_SIZE_FLOAT, vecVPluDvDt); plib_store_ps(outMatVecV + (i*nbVecCol + j)*PLIB_VECTOR_SIZE_FLOAT, vecVPlusDvDt);
// float uvSquare(u*v*v); // float uvSquare(u*v*v);
// float du(diffusionRateU*fullU - uvSquare + feedRate*(1.0f - u)); // float du(diffusionRateU*fullU - uvSquare + feedRate*(1.0f - u));
......
...@@ -7,4 +7,9 @@ set_property(TARGET naive_gray_scott PROPERTY COMPILE_FLAGS "-O3") ...@@ -7,4 +7,9 @@ set_property(TARGET naive_gray_scott PROPERTY COMPILE_FLAGS "-O3")
target_link_libraries(naive_gray_scott gray_scott_data_format gray_scott_naive tensor_alloc option_parser data_stream string_utils ${HDF5_CXX_LIBRARIES} TBB::tbb) target_link_libraries(naive_gray_scott gray_scott_data_format gray_scott_naive tensor_alloc option_parser data_stream string_utils ${HDF5_CXX_LIBRARIES} TBB::tbb)
add_executable(intrinsics_gray_scott main_intrinsics.cpp)
set_property(TARGET intrinsics_gray_scott PROPERTY COMPILE_FLAGS "-O3")
target_link_libraries(intrinsics_gray_scott gray_scott_data_format gray_scott_intrinsics tensor_alloc option_parser data_stream string_utils ${HDF5_CXX_LIBRARIES} TBB::tbb)
...@@ -82,12 +82,12 @@ bool simulateImage(size_t nbRow, size_t nbCol, size_t nbImage, size_t nbExtraSte ...@@ -82,12 +82,12 @@ bool simulateImage(size_t nbRow, size_t nbCol, size_t nbImage, size_t nbExtraSte
grayscott_propagation(tmpU2, tmpV2, tmpU1, tmpV1, nbRow, nbCol, grayscott_propagation(tmpU2, tmpV2, tmpU1, tmpV1, nbRow, nbCol,
matDeltaSquare, nbStencilRow, nbStencilCol, matDeltaSquare, nbStencilRow, nbStencilCol,
diffudionRateU, diffusionRateV, feedRate, killRate, dt); diffudionRateU, diffusionRateV, feedRate, killRate, dt);
fullMat.setRow(i, tmpV2);
///Let's swap the pointer ///Let's swap the pointer
swapValue(tmpU1, tmpU2); swapValue(tmpU1, tmpU2);
swapValue(tmpV1, tmpV2); swapValue(tmpV1, tmpV2);
} }
fullMat.setRow(i, tmpV2);
} }
progress.finish(); progress.finish();
std::cerr << "Done" << std::endl; std::cerr << "Done" << std::endl;
......
/***************************************
Auteur : Pierre Aubert
Mail : aubertp7@gmail.com
Licence : CeCILL-C
****************************************/
#include "OptionParser.h"
#include "temporary_alloc.h"
#include "ProgressTime.h"
#include "intrinsics_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("intrinsics_gray_scott --killrate=0.062 --feedrate=0.03 --nbimage=100 --nbrow=1080 --nbcol=1920 --output=outputFile.hdf5");
parser.setExampleShortOption("intrinsics_gray_scott -k 0.062 -f 0.03 -n 100 -r 1080 -c 1920 -o outputFile.hdf5");
float killRate(0.054f), feedRate(0.014f);
size_t nbImage(100lu), nbRow(100lu), nbCol(200lu);
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;
}
///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 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, 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);
fullMat.resize(nbImage);
PTensor<float> tmpInU, tmpInV, tmpOutU, tmpOutV;
float *tmpU1 = NULL, *tmpU2 = NULL, *tmpV1 = NULL, *tmpV2 = NULL;
allocate_temporary(tmpU1, tmpU2, tmpV1, tmpV2, tmpInU, tmpInV, tmpOutU, tmpOutV, nbRow, nbCol);
long nbStencilRow(3l), nbStencilCol(3l);
float diffudionRateU(0.1f), diffusionRateV(0.05f);
//This matrix of neigbour exchange is quite accurate but gives not so fun results
// float matDeltaSquare[] = {0.05f, 0.2f, 0.05f,
// 0.2f, 0.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};
//Let's convert these temporaries into intrinsics temporaries
PTensor<float> tmpVecInU(AllocMode::ALIGNED), tmpVecInV(AllocMode::ALIGNED), tmpVecOutU(AllocMode::ALIGNED), tmpVecOutV(AllocMode::ALIGNED);
tmpVecInU.fromScalToVecNeigbhour(tmpInU, PLIB_VECTOR_SIZE_FLOAT);
tmpVecInV.fromScalToVecNeigbhour(tmpInV, PLIB_VECTOR_SIZE_FLOAT);
tmpVecOutU.fromScalToVecNeigbhour(tmpOutU, PLIB_VECTOR_SIZE_FLOAT);
tmpVecOutV.fromScalToVecNeigbhour(tmpOutV, PLIB_VECTOR_SIZE_FLOAT);
PTensor<float> vecMatDeltaSquare(AllocMode::ALIGNED, nbStencilRow, nbStencilCol*PLIB_VECTOR_SIZE_FLOAT);
reshuffle_broadcastTensor(vecMatDeltaSquare.getData(), matDeltaSquare, nbStencilRow, nbStencilCol, 0lu, PLIB_VECTOR_SIZE_FLOAT);
tmpU1 = tmpVecInU.getData();
tmpU2 = tmpVecInV.getData();
tmpV1 = tmpVecOutU.getData();
tmpV2 = tmpVecOutV.getData();
float * ptrVecMatStencil = vecMatDeltaSquare.getData();
size_t nbVecRow(tmpVecInV.getFullNbRow()), nbVecCol(tmpVecInV.getNbCol());
PTensor<float> tmpScalOutV(AllocMode::ALIGNED);
ProgressTime progress(nbImage);
progress.start();
for(size_t i(0lu); i < nbImage; ++i){
progress.print();
for(size_t j(0lu); j < nbExtraStep; ++j){
grayscott_propagation(tmpU2, tmpV2, tmpU1, tmpV1, nbVecRow, nbVecCol, 0lu,
ptrVecMatStencil, nbStencilRow, nbStencilCol,
diffudionRateU, diffusionRateV, feedRate, killRate, dt);
//Let's update the dupplicated values
reshuffle_averageDupplicateVecNeighbour(tmpU2, nbVecRow, nbVecCol);
reshuffle_averageDupplicateVecNeighbour(tmpV2, nbVecRow, nbVecCol);
///Let's swap the pointer
swapValue(tmpU1, tmpU2);
swapValue(tmpV1, tmpV2);
}
if(tmpV2 == tmpVecOutV.getData()){
tmpScalOutV.fromVecToScalNeigbhour(tmpVecOutV);
}else{
tmpScalOutV.fromVecToScalNeigbhour(tmpVecInV); //The pointers were swaped
}
// fullMat.setRow(i, tmpScalOutV.getData());
// fullMat.setRow(i, tmpV1);
fullMat.setRow(i, tmpV2);
}
progress.finish();
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), 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, nbExtraStep, killRate, feedRate, dt, outputFile));
return b - 1;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment