From 883a0ddf2b12130b4e5daf54b3734d116b7e42a9 Mon Sep 17 00:00:00 2001 From: Yann Kempf Date: Tue, 17 Jan 2023 14:08:50 +0200 Subject: [PATCH 1/3] Added support for modern vg_ and fg_ variables in vlsvextract for rotation and plasma frame. --- tools/vlsvextract.cpp | 82 ++++++++++++++++++++++++++++++++----------- 1 file changed, 61 insertions(+), 21 deletions(-) diff --git a/tools/vlsvextract.cpp b/tools/vlsvextract.cpp index 18d8f2a07..a664b5db4 100644 --- a/tools/vlsvextract.cpp +++ b/tools/vlsvextract.cpp @@ -492,31 +492,44 @@ void getBulkVelocity(Real* V_bulk,vlsvinterface::Reader& vlsvReader,const string exit(1); } - // Read number density - double numberDensity; - double* ptr = &numberDensity; + // Read velocity + double velocity[3]; + double* ptr = velocity; xmlAttributes.clear(); xmlAttributes.push_back(make_pair("mesh",meshName)); - xmlAttributes.push_back(make_pair("name","rho")); + xmlAttributes.push_back(make_pair("name","proton/vg_v")); 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; - exit(1); - } + cerr << "Could not read velocity in " << __FILE__ << ":" << __LINE__ << ", trying to read density + momentum" << endl; + // Try old style + // 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; + exit(1); + } - V_bulk[0] = momentum[0] / (numberDensity + numeric_limits::min()); - V_bulk[1] = momentum[1] / (numberDensity + numeric_limits::min()); - V_bulk[2] = momentum[2] / (numberDensity + numeric_limits::min()); + V_bulk[0] = momentum[0] / (numberDensity + numeric_limits::min()); + V_bulk[1] = momentum[1] / (numberDensity + numeric_limits::min()); + V_bulk[2] = momentum[2] / (numberDensity + numeric_limits::min()); + } + V_bulk[0] = velocity[0]; + V_bulk[1] = velocity[1]; + V_bulk[2] = velocity[2]; } void getB(Real* B,vlsvinterface::Reader& vlsvReader,const string& meshName,const uint64_t& cellID) { @@ -570,6 +583,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 @@ -586,6 +600,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(); From a0d1d3ed2f6813701559bffe172d0144c8c7436d Mon Sep 17 00:00:00 2001 From: Yann Kempf Date: Tue, 17 Jan 2023 14:28:48 +0200 Subject: [PATCH 2/3] Multipop query of velocity for plasma frame in vlsvextract. --- tools/vlsvextract.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/vlsvextract.cpp b/tools/vlsvextract.cpp index a664b5db4..494c4976a 100644 --- a/tools/vlsvextract.cpp +++ b/tools/vlsvextract.cpp @@ -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; @@ -497,7 +497,7 @@ void getBulkVelocity(Real* V_bulk,vlsvinterface::Reader& vlsvReader,const string double* ptr = velocity; xmlAttributes.clear(); xmlAttributes.push_back(make_pair("mesh",meshName)); - xmlAttributes.push_back(make_pair("name","proton/vg_v")); + xmlAttributes.push_back(make_pair("name",popName+"/vg_v")); if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,ptr,false) == false) { cerr << "Could not read velocity in " << __FILE__ << ":" << __LINE__ << ", trying to read density + momentum" << endl; // Try old style @@ -736,7 +736,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); } From ce103c29e621f78913054f0887a3226a404a6c00 Mon Sep 17 00:00:00 2001 From: Yann Kempf Date: Tue, 17 Jan 2023 15:03:40 +0200 Subject: [PATCH 3/3] More complete set of trials to retrieve a velocity for plasma frame shifting. --- tools/vlsvextract.cpp | 83 ++++++++++++++++++++++++++++--------------- 1 file changed, 55 insertions(+), 28 deletions(-) diff --git a/tools/vlsvextract.cpp b/tools/vlsvextract.cpp index 494c4976a..801d73a22 100644 --- a/tools/vlsvextract.cpp +++ b/tools/vlsvextract.cpp @@ -492,44 +492,71 @@ void getBulkVelocity(Real* V_bulk,vlsvinterface::Reader& vlsvReader,const string exit(1); } - // Read velocity - double velocity[3]; - double* ptr = velocity; - 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) == false) { - cerr << "Could not read velocity in " << __FILE__ << ":" << __LINE__ << ", trying to read density + momentum" << 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 /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 /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; - double* ptr = &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) == false) { - cerr << "Could not read number density in " << __FILE__ << ":" << __LINE__ << endl; - exit(1); + 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::min()); + V_bulk[1] = momentum[1] / (numberDensity + numeric_limits::min()); + V_bulk[2] = momentum[2] / (numberDensity + numeric_limits::min()); + break; + } } - - // Read number density times velocity - double momentum[3]; - ptr = momentum; + // 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","rho_v")); - if (vlsvReader.read("VARIABLE",xmlAttributes,cellIndex,1,ptr,false) == false) { - cerr << "Could not read momentum in " << __FILE__ << ":" << __LINE__ << endl; - exit(1); + 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; } - - V_bulk[0] = momentum[0] / (numberDensity + numeric_limits::min()); - V_bulk[1] = momentum[1] / (numberDensity + numeric_limits::min()); - V_bulk[2] = momentum[2] / (numberDensity + numeric_limits::min()); - } - V_bulk[0] = velocity[0]; - V_bulk[1] = velocity[1]; - V_bulk[2] = velocity[2]; + // 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); + break; + } while (true); } void getB(Real* B,vlsvinterface::Reader& vlsvReader,const string& meshName,const uint64_t& cellID) {