/***************************************************************************************** * * * OpenSpace * * * * Copyright (c) 2014-2019 * * * * Permission is hereby granted, free of charge, to any person obtaining a copy of this * * software and associated documentation files (the "Software"), to deal in the Software * * without restriction, including without limitation the rights to use, copy, modify, * * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to * * permit persons to whom the Software is furnished to do so, subject to the following * * conditions: * * * * The above copyright notice and this permission notice shall be included in all copies * * or substantial portions of the Software. * * * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, * * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A * * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT * * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF * * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE * * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. * ****************************************************************************************/ #include #include #include #include #include #include #ifdef OPENSPACE_MODULE_KAMELEON_ENABLED #ifdef _MSC_VER #pragma warning (push) // Boost throws #pragma warning: there is no warning number '4675' #pragma warning (disable : 4619) #endif // _MSC_VER #include #include #include #ifdef _MSC_VER #pragma warning (pop) #endif // _MSC_VER #endif // OPENSPACE_MODULE_KAMELEON_ENABLED namespace { constexpr const char* _loggerCat = "FieldlinesSequence[ Kameleon ]"; constexpr const char* TAsPOverRho = "T = p/rho"; constexpr const char* JParallelB = "Current: mag(J||B)"; // [nPa]/[amu/cm^3] * ToKelvin => Temperature in Kelvin constexpr const float ToKelvin = 72429735.6984f; } // namespace namespace openspace::fls { // -------------------- DECLARE FUNCTIONS USED (ONLY) IN THIS FILE -------------------- // #ifdef OPENSPACE_MODULE_KAMELEON_ENABLED bool addLinesToState(ccmc::Kameleon* kameleon, const std::vector& seeds, const std::string& tracingVar, FieldlinesState& state); void addExtraQuantities(ccmc::Kameleon* kameleon, std::vector& extraScalarVars, std::vector& extraMagVars, FieldlinesState& state); void prepareStateAndKameleonForExtras(ccmc::Kameleon* kameleon, std::vector& extraScalarVars, std::vector& extraMagVars, FieldlinesState& state); #endif // OPENSPACE_MODULE_KAMELEON_ENABLED // ------------------------------------------------------------------------------------ // /** Traces field lines from the provided cdf file using kameleon and stores the data in * the provided FieldlinesState. * Returns `false` if it fails to create a valid state. Requires the kameleon module to * be activated! * \param state, FieldlineState which should hold the extracted data * \param cdfPath, std::string of the absolute path to a .cdf file * \param seedPoints, vector of seed points from which to trace field lines * \param tracingVar, which quantity to trace lines from. Typically "b" for magnetic field * lines and "u" for velocity flow lines * \param extraVars, extra scalar quantities to be stored in the FieldlinesState; e.g. "T" * for temperature, "rho" for density or "P" for pressure * \param extraMagVars, variables which should be used for extracting magnitudes, must be * a multiple of 3; e.g. "ux", "uy" & "uz" to get the magnitude of the velocity * vector at each line vertex */ bool convertCdfToFieldlinesState(FieldlinesState& state, const std::string& cdfPath, const std::vector& seedPoints, const std::string& tracingVar, std::vector& extraVars, std::vector& extraMagVars) { #ifndef OPENSPACE_MODULE_KAMELEON_ENABLED LERROR("CDF inputs provided but Kameleon module is deactivated"); return false; #else // OPENSPACE_MODULE_KAMELEON_ENABLED // Create Kameleon object and open CDF file! std::unique_ptr kameleon = kameleonHelper::createKameleonObject( cdfPath ); state.setModel(fls::stringToModel(kameleon->getModelName())); state.setTriggerTime(kameleonHelper::getTime(kameleon.get())); if (addLinesToState(kameleon.get(), seedPoints, tracingVar, state)) { // The line points are in their RAW format (unscaled & maybe spherical) // Before we scale to meters (and maybe cartesian) we must extract // the extraQuantites, as the iterpolator needs the unaltered positions addExtraQuantities(kameleon.get(), extraVars, extraMagVars, state); switch (state.model()) { case fls::Model::Batsrus: state.scalePositions(fls::ReToMeter); break; case fls::Model::Enlil: state.convertLatLonToCartesian(fls::AuToMeter); break; default: break; } return true; } return false; #endif // OPENSPACE_MODULE_KAMELEON_ENABLED } #ifdef OPENSPACE_MODULE_KAMELEON_ENABLED /** * Traces and adds line vertices to state. * Vertices are not scaled to meters nor converted from spherical into cartesian * coordinates. * Note that extraQuantities will NOT be set! */ bool addLinesToState(ccmc::Kameleon* kameleon, const std::vector& seedPoints, const std::string& tracingVar, FieldlinesState& state) { float innerBoundaryLimit; switch (state.model()) { case fls::Model::Batsrus: innerBoundaryLimit = 2.5f; // TODO specify in Lua? break; case fls::Model::Enlil: innerBoundaryLimit = 0.11f; // TODO specify in Lua? break; default: LERROR( "OpenSpace's fieldlines sequence currently only supports CDFs from the " "BATSRUS and ENLIL models!" ); return false; } // ---------------------------- LOAD TRACING VARIABLE ---------------------------- // if (!kameleon->loadVariable(tracingVar)) { LERROR("Failed to load tracing variable: " + tracingVar); return false; } bool success = false; LINFO("Tracing field lines!"); // LOOP THROUGH THE SEED POINTS, TRACE LINES, CONVERT POINTS TO glm::vec3 AND STORE // for (const glm::vec3& seed : seedPoints) { //--------------------------------------------------------------------------// // We have to create a new tracer (or actually a new interpolator) for each // // new line, otherwise some issues occur // //--------------------------------------------------------------------------// std::unique_ptr interpolator = std::make_unique(kameleon->model); ccmc::Tracer tracer(kameleon, interpolator.get()); tracer.setInnerBoundary(innerBoundaryLimit); // TODO specify in Lua? ccmc::Fieldline ccmcFieldline = tracer.bidirectionalTrace( tracingVar, seed.x, seed.y, seed.z ); const std::vector& positions = ccmcFieldline.getPositions(); const size_t nLinePoints = positions.size(); std::vector vertices; vertices.reserve(nLinePoints); for (const ccmc::Point3f& p : positions) { vertices.emplace_back(p.component1, p.component2, p.component3); } state.addLine(vertices); success |= (nLinePoints > 0); } return success; } #endif // OPENSPACE_MODULE_KAMELEON_ENABLED /** * Loops through state's _vertexPositions and extracts corresponding 'extraQuantities' * from the kameleon object using a ccmc::interpolator. * Note that the positions MUST be unaltered (NOT scaled NOR converted to a different * coordinate system)! * * @param kameleon raw pointer to an already opened Kameleon object * @param extraScalarVars vector of strings. Strings should be names of a scalar * quantities to load into _extraQuantites; such as: "T" for temperature or "rho" for * density. * @param extraMagVars vector of strings. Size must be multiple of 3. Strings should be * names of the components needed to calculate magnitude. E.g. {"ux", "uy", "uz"} will * calculate: sqrt(ux*ux + uy*uy + uz*uz). Magnitude will be stored in _extraQuantities * @param state, The FieldlinesState which the extra quantities should be added to. */ #ifdef OPENSPACE_MODULE_KAMELEON_ENABLED void addExtraQuantities(ccmc::Kameleon* kameleon, std::vector& extraScalarVars, std::vector& extraMagVars, FieldlinesState& state) { prepareStateAndKameleonForExtras(kameleon, extraScalarVars, extraMagVars, state); const size_t nXtraScalars = extraScalarVars.size(); const size_t nXtraMagnitudes = extraMagVars.size() / 3; std::unique_ptr interpolator = std::make_unique(kameleon->model); // ------ Extract all the extraQuantities from kameleon and store in state! ------ // for (const glm::vec3& p : state.vertexPositions()) { // Load the scalars! for (size_t i = 0; i < nXtraScalars; i++) { float val; if (extraScalarVars[i] == TAsPOverRho) { val = interpolator->interpolate("p", p.x, p.y, p.z); val *= ToKelvin; val /= interpolator->interpolate("rho", p.x, p.y, p.z); } else { val = interpolator->interpolate(extraScalarVars[i], p.x, p.y, p.z); // When measuring density in ENLIL CCMC multiply by the radius^2 if (extraScalarVars[i] == "rho" && state.model() == fls::Model::Enlil) { val *= std::pow(p.x * fls::AuToMeter, 2.0f); } } state.appendToExtra(i, val); } // Calculate and store the magnitudes! for (size_t i = 0; i < nXtraMagnitudes; ++i) { const size_t idx = i*3; const float x = interpolator->interpolate(extraMagVars[idx] , p.x, p.y, p.z); const float y = interpolator->interpolate(extraMagVars[idx+1], p.x, p.y, p.z); const float z = interpolator->interpolate(extraMagVars[idx+2], p.x, p.y, p.z); float val; // When looking at the current's magnitude in Batsrus, CCMC staff are // only interested in the magnitude parallel to the magnetic field if (state.extraQuantityNames()[nXtraScalars + i] == JParallelB) { const glm::vec3 normMagnetic = glm::normalize(glm::vec3( interpolator->interpolate("bx", p.x, p.y, p.z), interpolator->interpolate("by", p.x, p.y, p.z), interpolator->interpolate("bz", p.x, p.y, p.z))); // Magnitude of the part of the current vector that's parallel to // the magnetic field vector! val = glm::dot(glm::vec3(x,y,z), normMagnetic); } else { val = std::sqrt(x*x + y*y + z*z); } state.appendToExtra(i + nXtraScalars, val); } } } #endif // OPENSPACE_MODULE_KAMELEON_ENABLED /** Validate the provided extra quantity variables -> load the data from the validated * quantities into the kameleon object & add the quantity names into the state's * _extraQuantityNames vector. * * @param kameleon, raw pointer to an already opened kameleon object * @param extraScalarVars, names of scalar quantities to add to state; e.g "rho" for * density * @param extraMagVars, names of the variables used for calculating magnitudes. Must be * multiple of 3. */ #ifdef OPENSPACE_MODULE_KAMELEON_ENABLED void prepareStateAndKameleonForExtras(ccmc::Kameleon* kameleon, std::vector& extraScalarVars, std::vector& extraMagVars, FieldlinesState& state) { std::vector extraQuantityNames; fls::Model model = fls::stringToModel(kameleon->getModelName()); // Load the existing SCALAR variables into kameleon. // Remove non-existing variables from vector for (int i = 0; i < static_cast(extraScalarVars.size()); ++i) { std::string& str = extraScalarVars[i]; bool success = kameleon->doesVariableExist(str) && kameleon->loadVariable(str); if (!success && (model == fls::Model::Batsrus && (str == TAsPOverRho || str == "T"))) { LDEBUG("BATSRUS doesn't contain variable T for temperature. Trying to " "calculate it using the ideal gas law: T = pressure/density"); constexpr const char* p = "p"; constexpr const char* r = "rho"; success = kameleon->doesVariableExist(p) && kameleon->loadVariable(p) && kameleon->doesVariableExist(r) && kameleon->loadVariable(r); str = TAsPOverRho; } if (!success) { LWARNING(fmt::format( "Failed to load extra variable: '{}'. Ignoring", str )); extraScalarVars.erase(extraScalarVars.begin() + i); --i; } else { extraQuantityNames.push_back(str); } } // Load the existing magnitude variables (should be provided in multiple of 3) // into kameleon. Remove non-existing variables from vector if (extraMagVars.size() % 3 == 0) { for (size_t i = 0; i < extraMagVars.size(); i += 3) { const std::string& s1 = extraMagVars[i]; const std::string& s2 = extraMagVars[i+1]; const std::string& s3 = extraMagVars[i+2]; bool success = kameleon->doesVariableExist(s1) && kameleon->doesVariableExist(s2) && kameleon->doesVariableExist(s3) && kameleon->loadVariable(s1) && kameleon->loadVariable(s2) && kameleon->loadVariable(s3); std::string name = "Magnitude of (" + s1 + ", " + s2 + ", " + s3 + ")"; if (success && model == fls::Model::Batsrus && s1 == "jx" && s2 == "jy" && s3 == "jz") { // CCMC isn't really interested in the magnitude of current, but by the // magnitude of the part of the current's vector that is parallel to the // magnetic field => ensure that the magnetic variables are loaded success = kameleon->doesVariableExist("bx") && kameleon->doesVariableExist("by") && kameleon->doesVariableExist("bz") && kameleon->loadVariable("bx") && kameleon->loadVariable("by") && kameleon->loadVariable("bz"); name = JParallelB; } if (!success) { LWARNING(fmt::format( "Failed to load at least one of the magnitude variables: {}, {}, {} " "& {}. Removing ability to store corresponding magnitude", s1, s2, s3 )); extraMagVars.erase( extraMagVars.begin() + i, extraMagVars.begin() + i + 3 ); i -= 3; } else { extraQuantityNames.push_back(std::move(name)); } } } else { // WRONG NUMBER OF MAGNITUDE VARIABLES.. REMOVE ALL! extraMagVars.clear(); LWARNING(fmt::format( "Wrong number of variables provided for storing magnitudes. Expects multiple " "of 3 but {} are provided", extraMagVars.size() )); } state.setExtraQuantityNames(std::move(extraQuantityNames)); } #endif // OPENSPACE_MODULE_KAMELEON_ENABLED } // namespace openspace::fls