Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Ionosphere vlsvdiff #2

Open
wants to merge 11 commits into
base: dev
Choose a base branch
from
Open

Ionosphere vlsvdiff #2

Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
daa99c9
vlsvdiff now works for AMR runs
kstppd May 12, 2021
03db7e7
Merge branch 'dev' of github.com:kstppd/vlasiator into dev
ykempf May 18, 2021
ea6f26e
vlsvdiff now works for AMR runs
kstppd May 18, 2021
b5b3d54
Merge branch 'dev' of github.com:kstppd/vlasiator into dev
ykempf May 18, 2021
5d1a139
increment by variable vector size
kstppd May 18, 2021
8b5f93e
Merge branch 'dev' of github.com:kstppd/vlasiator into dev
ykempf May 18, 2021
f6b92ac
remove some print that was left in there by accident
kstppd May 19, 2021
190e833
Merge branch 'dev' of github.com:kstppd/vlasiator into dev
ykempf May 19, 2021
2d833ba
This will handle datasize mismatches
kstppd May 20, 2021
ef5733f
Merge branch 'dev' of github.com:kstppd/vlasiator into ionosphere_vls…
ykempf May 20, 2021
b3a852c
vlsvdiff for ionosphere
ykempf May 20, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
204 changes: 135 additions & 69 deletions tools/vlsvdiff.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -421,14 +421,13 @@ bool convertMesh(vlsvinterface::Reader& vlsvReader,
cerr << "ERROR, failed to get array info for '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
return false;
}

if (gridName==gridType::SpatialGrid) {
char *variableBuffer = new char[variableVectorSize * variableDataSize];
float *variablePtrFloat = reinterpret_cast<float *>(variableBuffer);
double *variablePtrDouble = reinterpret_cast<double *>(variableBuffer);
uint *variablePtrUint = reinterpret_cast<uint *>(variableBuffer);
int *variablePtrInt = reinterpret_cast<int *>(variableBuffer);


if (gridName==gridType::SpatialGrid){

// Read the mesh array one node (of a spatial cell) at a time
// and create a map which contains each cell's CellID and variable to be extracted
Expand Down Expand Up @@ -488,32 +487,49 @@ bool convertMesh(vlsvinterface::Reader& vlsvReader,
if (storeCellOrder == true) {
cellOrder[CellID] = i;
}
}

}else if (gridName==gridType::fsgrid){


}
delete variableBuffer;
} else if (gridName==gridType::fsgrid) {


// Get Spatial Grid's max refinement Level
int maxRefLevel=0;
list<pair<string, string>> meshAttributesIn;
meshAttributesIn.push_back(make_pair("name", "SpatialGrid"));
map<string,string> meshAttributesOut;
if (vlsvReader.getArrayAttributes("MESH", meshAttributesIn,meshAttributesOut) == false)
{
cerr << "ERROR, failed to get array info for '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
return false;
}

std::map<string, string>::iterator attributesOutIt;
attributesOutIt = meshAttributesOut.find("max_refinement_level");
if (attributesOutIt != meshAttributesOut.end())
{
maxRefLevel = stoi(attributesOutIt->second);
}
int numtasks;
int xcells,ycells,zcells;
int xcells,ycells,zcells;
vlsvReader.readParameter("numWritingRanks",numtasks);
vlsvReader.readParameter("xcells_ini",xcells);
vlsvReader.readParameter("ycells_ini",ycells);
vlsvReader.readParameter("zcells_ini",zcells);
xcells*=pow(2,maxRefLevel);
ycells*=pow(2,maxRefLevel);
zcells*=pow(2,maxRefLevel);
std::array<int,3> GlobalBox={xcells,ycells,zcells};
std::array<int,3> thisDomainDecomp;

//Compute Domain Decomposition Scheme for this vlsv file
computeDomainDecomposition(GlobalBox,numtasks,thisDomainDecomp);


std::array<int32_t,3> taskSize,taskStart;
std::array<int32_t,3> taskEnd;
int readOffset=0;
size_t readSize;
int index,my_x,my_y,my_z;
orderedData->clear();

//Read into buffer
for (int task=0; task<numtasks; task++){

my_x=task/thisDomainDecomp[2]/thisDomainDecomp[1];
Expand All @@ -533,54 +549,109 @@ bool convertMesh(vlsvinterface::Reader& vlsvReader,
taskEnd[1]= taskStart[1]+taskSize[1];
taskEnd[2]= taskStart[2]+taskSize[2];

readSize= taskSize[0] * taskSize[1] * taskSize[2];
std::vector<Real> readIn(variableVectorSize * variableDataSize*readSize);


int counter2=0;
uint64_t globalindex;
int64_t counter=0;
for(int z=taskStart[2]; z<taskEnd[2]; z++) {
for(int y=taskStart[1]; y<taskEnd[1]; y++) {
for(int x=taskStart[0]; x<taskEnd[0]; x++) {

//Get global index
globalindex= x + y*xcells + z*xcells*ycells;

if (vlsvReader.readArray("VARIABLE", variableAttributes, readOffset+counter,1, variableBuffer) == false) {
cerr << "ERROR, failed to read variable '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
variableSuccess = false;
abort();
break;
}

// Get the variable value
Real extract = NAN;
int64_t readSize= taskSize[0] * taskSize[1] * taskSize[2] ;
//Allocate vector for reading
std::vector<Real> buffer(readSize*variableVectorSize);

if ( variableDataSize==sizeof(Real)){
if (vlsvReader.readArray("VARIABLE", variableAttributes, readOffset, readSize, (char*)buffer.data()) == false) {
cerr << "ERROR, failed to read variable '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
variableSuccess = false;
break;
}
}else{
std::vector<float> tmpbuffer(readSize * variableVectorSize);
if (vlsvReader.readArray("VARIABLE", variableAttributes, readOffset, readSize, (char *)tmpbuffer.data()) == false){
cerr << "ERROR, failed to read variable '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
variableSuccess = false;
break;
}
for (int i = 0; i < readSize * variableVectorSize; i++){
buffer[i] = tmpbuffer[i];
}
}

switch (variableDataType) {
uint64_t globalindex,counter=0;;
for (int z=taskStart[2]; z<taskEnd[2]; z++){
for (int y=taskStart[1]; y< taskEnd[1]; y++){
for (int x=taskStart[0]; x<taskEnd[0]; x++){
globalindex= x + y*xcells + z*xcells*ycells;
Real data;
switch (variableDataType){
case datatype::type::FLOAT:
if(variableDataSize == sizeof(float)) extract = (Real)(variablePtrFloat[compToExtract]);
if(variableDataSize == sizeof(double)) extract = (Real)(variablePtrDouble[compToExtract]);
if (variableDataSize == sizeof(float))
memcpy(&data, &buffer[counter + compToExtract], sizeof(float));
if (variableDataSize == sizeof(double))
memcpy(&data, &buffer[counter + compToExtract], sizeof(double));
break;
case datatype::type::UINT:
extract = (Real)(variablePtrUint[compToExtract]);
memcpy(&data, &buffer[counter + compToExtract], sizeof(uint));
break;
case datatype::type::INT:
extract = (Real)(variablePtrInt[compToExtract]);
memcpy(&data, &buffer[counter + compToExtract], sizeof(int));
break;
case datatype::type::UNKNOWN:
cerr << "ERROR, BAD DATATYPE AT " << __FILE__ << " " << __LINE__ << endl;
break;
}
//Add to map
orderedData->insert(pair<uint64_t, Real>(globalindex, data));
counter+=variableVectorSize;
}
orderedData->insert(pair<uint64_t, Real>(globalindex, extract));
counter++;

}
}
}
readOffset+=readSize;
}
}else{
} else if (gridName==gridType::ionosphere) {
char *variableBuffer = new char[variableArraySize * variableVectorSize * variableDataSize];
float *variablePtrFloat = reinterpret_cast<float *>(variableBuffer);
double *variablePtrDouble = reinterpret_cast<double *>(variableBuffer);
uint *variablePtrUint = reinterpret_cast<uint *>(variableBuffer);
int *variablePtrInt = reinterpret_cast<int *>(variableBuffer);

if (compToExtract + 1 > variableVectorSize) {
cerr << "ERROR invalid component, this variable has size " << variableVectorSize << endl;
abort();
}

if (storeCellOrder == true) {
cellOrder.clear();
}

orderedData->clear();

if (vlsvReader.readArray("VARIABLE", variableAttributes, 0, variableArraySize, variableBuffer) == false) {
cerr << "ERROR, failed to read variable '" << _varToExtract << "' at " << __FILE__ << " " << __LINE__ << endl;
variableSuccess = false;
}

for (uint64_t i=0; i<variableArraySize; ++i) {
// Get the variable value
Real extract = NAN;

switch (variableDataType) {
case datatype::type::FLOAT:
if(variableDataSize == sizeof(float)) extract = (Real)(variablePtrFloat[i*variableVectorSize + compToExtract]);
if(variableDataSize == sizeof(double)) extract = (Real)(variablePtrDouble[i*variableVectorSize + compToExtract]);
break;
case datatype::type::UINT:
extract = (Real)(variablePtrUint[i*variableVectorSize + compToExtract]);
break;
case datatype::type::INT:
extract = (Real)(variablePtrInt[i*variableVectorSize + compToExtract]);
break;
case datatype::type::UNKNOWN:
cerr << "ERROR, BAD DATATYPE AT " << __FILE__ << " " << __LINE__ << endl;
break;
}
// Put those into the map
orderedData->insert(pair<uint64_t, Real>(i, extract));
if (storeCellOrder == true) {
cellOrder[i] = i;
}
}
delete variableBuffer;
} else {
cerr<<"meshName not recognized\t" << __FILE__ << " " << __LINE__ <<endl;
abort();
}
Expand All @@ -591,7 +662,6 @@ bool convertMesh(vlsvinterface::Reader& vlsvReader,
if (variableSuccess == false) {
cerr << "ERROR reading array VARIABLE " << varToExtract << endl;
}
delete variableBuffer;
return meshSuccess && variableSuccess;
}

Expand Down Expand Up @@ -737,14 +807,12 @@ bool pDistance(const map<uint, Real>& orderedData1,
value = abs(it1->second - it2->second);
*absolute = max(*absolute, value);
length = max(length, abs(it1->second));

}
if (gridName==gridType::SpatialGrid){
}
if (gridName==gridType::SpatialGrid) {
array[cellOrder.at(it1->first)] = value;
}else if (gridName==gridType::fsgrid) {

} else if (gridName==gridType::fsgrid || gridName==gridType::ionosphere) {
array.at(it1->first)=value;
}
}
}
} else if (p == 1) {
for (map<uint,Real>::const_iterator it1=orderedData1.begin(); it1!=orderedData1.end(); ++it1) {
Expand All @@ -754,13 +822,12 @@ bool pDistance(const map<uint, Real>& orderedData1,
value = abs(it1->second - it2->second);
*absolute += value;
length += abs(it1->second);

}
if (gridName==gridType::SpatialGrid){
}
if (gridName==gridType::SpatialGrid) {
array[cellOrder.at(it1->first)] = value;
}else if (gridName==gridType::fsgrid){
} else if (gridName==gridType::fsgrid || gridName==gridType::ionosphere) {
array[it1->first]=value;
}
}
}
} else {
for (map<uint,Real>::const_iterator it1=orderedData1.begin(); it1!=orderedData1.end(); ++it1) {
Expand All @@ -770,13 +837,12 @@ bool pDistance(const map<uint, Real>& orderedData1,
value = pow(abs(it1->second - it2->second), p);
*absolute += value;
length += pow(abs(it1->second), p);

}
if (gridName==gridType::SpatialGrid){
}
if (gridName==gridType::SpatialGrid) {
array[cellOrder.at(it1->first)] = pow(value,1.0/p);
}else if (gridName==gridType::fsgrid){
}else if (gridName==gridType::fsgrid || gridName==gridType::ionosphere) {
array[it1->first]=pow(value,1.0/p);
}
}
}
*absolute = pow(*absolute, 1.0 / p);
length = pow(length, 1.0 / p);
Expand Down Expand Up @@ -1788,13 +1854,13 @@ int main(int argn,char* args[]) {


//Figure out Meshname
if (attributes["--meshname"] == "SpatialGrid") {
gridName=gridType::SpatialGrid ;
}else if (attributes["--meshname"]=="fsgrid"){
gridName=gridType::fsgrid ;
}else if (attributes["--meshname"]=="ionosphere"){
gridName=gridType::ionosphere ;
}else{
if (attributes["--meshname"] == "SpatialGrid") {
gridName=gridType::SpatialGrid;
} else if (attributes["--meshname"]=="fsgrid") {
gridName=gridType::fsgrid;
} else if (attributes["--meshname"]=="ionosphere") {
gridName=gridType::ionosphere;
} else {
std::cout<<attributes["--meshname"]<<std::endl;
std::cerr<<"Wrong grid type"<<std::endl;
abort();
Expand Down