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

Added support for modern vg_ and fg_ variables in vlsvextract for rotation and plasma frame #689

Merged
merged 3 commits into from
Feb 9, 2023
Merged

Added support for modern vg_ and fg_ variables in vlsvextract for rotation and plasma frame #689

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
883a0dd
Added support for modern vg_ and fg_ variables in vlsvextract for rot…
ykempf Jan 17, 2023
a0d1d3e
Multipop query of velocity for plasma frame in vlsvextract.
ykempf Jan 17, 2023
ce103c2
More complete set of trials to retrieve a velocity for plasma frame s…
ykempf Jan 17, 2023
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
119 changes: 93 additions & 26 deletions tools/vlsvextract.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -447,7 +447,7 @@ void applyRotation(const Real* B,Real* transform) {
}
}

void getBulkVelocity(Real* V_bulk,vlsvinterface::Reader& vlsvReader,const string& meshName,const uint64_t& cellID) {
void getBulkVelocity(Real* V_bulk,vlsvinterface::Reader& vlsvReader,const string& meshName,const string& popName,const uint64_t& cellID) {
//Declarations
vlsv::datatype::type cellIdDataType;
uint64_t cellIdArraySize, cellIdVectorSize, cellIdDataSize;
Expand Down Expand Up @@ -492,31 +492,71 @@ void getBulkVelocity(Real* V_bulk,vlsvinterface::Reader& vlsvReader,const string
exit(1);
}

// Read number density
double numberDensity;
double* ptr = &numberDensity;
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name","rho"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,ptr,false) == false) {
cerr << "Could not read number density in " << __FILE__ << ":" << __LINE__ << endl;
exit(1);
}

// Read number density times velocity
double momentum[3];
ptr = momentum;
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name","rho_v"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,ptr,false) == false) {
cerr << "Could not read momentum in " << __FILE__ << ":" << __LINE__ << endl;
do {
// Read combined vg_v
double velocity[3];
double* ptr = velocity;
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name","vg_v"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,ptr,false) == true) {
cerr << "NOTE: Using combined vg_v (all populations!) for plasma frame shifting." << endl;
V_bulk[0] = velocity[0];
V_bulk[0] = velocity[0];
V_bulk[0] = velocity[0];
break;
}
// Read <pop>/vg_v
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name",popName+"/vg_v"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,ptr,false) == true) {
cerr << "NOTE: Using <pop>/vg_v for plasma frame shifting." << endl;
V_bulk[0] = velocity[0];
V_bulk[0] = velocity[0];
V_bulk[0] = velocity[0];
break;
}
// Try old style
// Read number density
double numberDensity;
ptr = &numberDensity;
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name","rho"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,ptr,false) == true) {
// Read number density times velocity
double momentum[3];
ptr = momentum;
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name","rho_v"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,ptr,false) == true) {
cerr << "NOTE: Using rho_v / rho for plasma frame shifting." << endl;
V_bulk[0] = momentum[0] / (numberDensity + numeric_limits<double>::min());
V_bulk[1] = momentum[1] / (numberDensity + numeric_limits<double>::min());
V_bulk[2] = momentum[2] / (numberDensity + numeric_limits<double>::min());
break;
}
}
// Read combined vg_v in restart file style
double moments[5];
ptr = moments;
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name","moments"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,ptr,false) == true) {
cerr << "NOTE: Using combined vg_v (all populations!) from restart for plasma frame shifting." << endl;
V_bulk[0] = moments[1];
V_bulk[0] = moments[2];
V_bulk[0] = moments[3];
break;
}
// We should have broken out before or we're doomed
cerr << "ERROR: Could not find a usable velocity for plasma frame shift!" << endl;
exit(1);
}

V_bulk[0] = momentum[0] / (numberDensity + numeric_limits<double>::min());
V_bulk[1] = momentum[1] / (numberDensity + numeric_limits<double>::min());
V_bulk[2] = momentum[2] / (numberDensity + numeric_limits<double>::min());
break;
} while (true);
}

void getB(Real* B,vlsvinterface::Reader& vlsvReader,const string& meshName,const uint64_t& cellID) {
Expand Down Expand Up @@ -570,6 +610,7 @@ void getB(Real* B,vlsvinterface::Reader& vlsvReader,const string& meshName,const

// Magnetic field can exists in the file in few different variables.
// Here we go with the following priority:
// - vg_b_vol
// - B_vol
// - BGB_vol + PERB_vol
// - B
Expand All @@ -586,6 +627,32 @@ void getB(Real* B,vlsvinterface::Reader& vlsvReader,const string& meshName,const

bool B_read = true;
do {
// Attempt to read 'vg_b_vol'
B_read = true;
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name","vg_b_vol"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,B1_ptr,false) == false) B_read = false;
if (B_read == true) {
if (runDebug == true) cerr << "Using vg_b_vol" << endl;
break;
}

// Attempt to read 'vg_b_background_vol' + 'vg_b_perturbed_vol'
B_read = true;
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name","vg_b_background_vol"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,B1_ptr,false) == false) B_read = false;
xmlAttributes.clear();
xmlAttributes.push_back(make_pair("mesh",meshName));
xmlAttributes.push_back(make_pair("name","vg_b_perturbed_vol"));
if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,B2_ptr,false) == false) B_read = false;
if (B_read == true) {
if (runDebug == true) cerr << "Using vg_b_background_vol + vg_b_perturbed_vol" << endl;
break;
}

// Attempt to read 'B_vol'
B_read = true;
xmlAttributes.clear();
Expand Down Expand Up @@ -696,7 +763,7 @@ bool convertVelocityBlocks2(

if (plasmaFrame == true) {
Real V_bulk[3];
getBulkVelocity(V_bulk,vlsvReader,meshName,cellID);
getBulkVelocity(V_bulk,vlsvReader,meshName,popName,cellID);
applyTranslation(V_bulk,transform);
}

Expand Down