Access solution and primitives via python wrapper by pcarruscag · Pull Request #1938 · su2code/SU2 · GitHub
Skip to content

Access solution and primitives via python wrapper #1938

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

Merged
bigfooted merged 3 commits into from
Mar 1, 2023
Merged

Access solution and primitives via python wrapper #1938

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
62aaee7
matrix view type to manipulate fields in python without copying data
pcarruscag Feb 26, 2023
e444b31
Merge remote-tracking branch 'upstream/develop' into python_access_so…
pcarruscag Feb 26, 2023
e296613
extra line and update docker images
pcarruscag Feb 26, 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
18 changes: 9 additions & 9 deletions .github/workflows/regression.yml
View file
Open in desktop
2 changes: 1 addition & 1 deletion .github/workflows/release-management.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ jobs:
key: ${{ matrix.os_bin }}-${{ github.sha }}
restore-keys: ${{ matrix.os_bin }}
- name: Build
uses: docker://ghcr.io/su2code/su2/build-su2-cross:221224-1158
uses: docker://ghcr.io/su2code/su2/build-su2-cross:230225-2136
with:
args: -b ${{ github.sha }} -f "${{matrix.flags}}"
- name: Create Archive
Expand Down
115 changes: 115 additions & 0 deletions Common/include/containers/CPyWrapperMatrixView.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,115 @@
/*!
* \file CPyWrapperMatrixView.hpp
* \brief A simple matrix view to use with the python wrapper.
* \author P. Gomes
* \version 7.5.1 "Blackbird"
*
* SU2 Project Website: https://su2code.github.io
*
* The SU2 Project is maintained by the SU2 Foundation
* (http://su2foundation.org)
*
* Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md)
*
* SU2 is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* SU2 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with SU2. If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once

#include <utility>
#include <vector>
#include <string>

#include "C2DContainer.hpp"
#include "../parallelization/mpi_structure.hpp"

/*!
* \class CPyWrapperMatrixView
* \ingroup PySU2
* \brief A simple matrix view to use with the python wrapper. The accessors
* in this class provide a passive double interface to su2activematrix.
* This can be extended to allow access to derivative information of su2double.
*/
class CPyWrapperMatrixView {
private:
static_assert(su2activematrix::IsRowMajor, "");
su2double* data_ = nullptr;
unsigned long rows_ = 0, cols_ = 0;
std::string name_;
bool read_only_ = false;

inline const su2double& Access(unsigned long row, unsigned long col) const {
if (row > rows_ || col > cols_) SU2_MPI::Error(name_ + " out of bounds", "CPyWrapperMatrixView");
return data_[row * cols_ + col];
}
inline su2double& Access(unsigned long row, unsigned long col) {
if (read_only_) SU2_MPI::Error(name_ + " is read-only", "CPyWrapperMatrixView");
const auto& const_me = *this;
return const_cast<su2double&>(const_me.Access(row, col));
}

public:
CPyWrapperMatrixView() = default;

/*!
* \brief Construct the view of the matrix.
* \note "name" should be set to the variable name being returned to give better information to users.
* \note "read_only" can be set to true to prevent the data from being modified.
*/
CPyWrapperMatrixView(su2activematrix& mat, const std::string& name, bool read_only)
: data_(mat.data()), rows_(mat.rows()), cols_(mat.cols()), name_(name), read_only_(read_only) {}

/*!
* \brief Returns the shape of the matrix.
*/
std::pair<unsigned long, unsigned long> Shape() const { return std::make_pair(rows_, cols_); }

/*!
* \brief Returns whether the data is read-only [true] or if it can be modified [false].
*/
bool IsReadOnly() const { return read_only_; }

/*!
* \brief Gets the value for a (row, column) pair.
*/
passivedouble operator() (unsigned long row, unsigned long col) const { return Get(row, col); }

/*!
* \brief Gets the value for a (row, column) pair.
*/
passivedouble Get(unsigned long row, unsigned long col) const { return SU2_TYPE::GetValue(Access(row, col)); }

/*!
* \brief Gets the values for a row of the matrix.
*/
std::vector<passivedouble> Get(unsigned long row) const {
std::vector<passivedouble> vals(cols_);
for (unsigned long j = 0; j < cols_; ++j) vals[j] = Get(row, j);
return vals;
}

/*!
* \brief Sets the value for a (row, column) pair.
* \note This clears derivative information (consistently with C++ operator= with passive rhs).
*/
void Set(unsigned long row, unsigned long col, passivedouble val) { Access(row, col) = val; }

/*!
* \brief Sets the values for a row of the matrix.
*/
void Set(unsigned long row, std::vector<passivedouble> vals) {
unsigned long j = 0;
for (const auto& val : vals) Set(row, j++, val);
}
};
2 changes: 1 addition & 1 deletion Common/src/containers/meson.build
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
common_src += files(['CTrapezoidalMap.cpp',
'CFileReaderLUT.cpp',
'CFileReaderLUT.cpp',
'CLookUpTable.cpp'])
64 changes: 48 additions & 16 deletions SU2_CFD/include/drivers/CDriverBase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
#pragma once

#include "../../../Common/include/CConfig.hpp"
#include "../../../Common/include/containers/CPyWrapperMatrixView.hpp"
#include "../numerics/CNumerics.hpp"
#include "../output/COutput.hpp"
#include "../solvers/CSolver.hpp"
Expand All @@ -38,7 +39,6 @@
* \brief Base class for all drivers.
* \author H. Patel, A. Gastaldi
*/

class CDriverBase {
protected:
int rank, /*!< \brief MPI Rank. */
Expand Down Expand Up @@ -189,25 +189,24 @@ class CDriverBase {
bool GetNodeDomain(unsigned long iPoint) const;

/*!
* \brief Get the initial (un-deformed) coordinates of a mesh node.
* \param[in] iPoint - Mesh node index.
* \return Initial node coordinates (nDim).
*/
vector<passivedouble> GetInitialCoordinates(unsigned long iPoint) const;

/*!
* \brief Get the current coordinates of a mesh node.
* \param[in] iPoint - Mesh node index.
* \return Node coordinates (nDim).
* \brief Get a read-only view of the initial (undeformed) coordinates of all mesh nodes.
*/
vector<passivedouble> GetCoordinates(unsigned long iPoint) const;
inline CPyWrapperMatrixView InitialCoordinates() const {
if (!main_config->GetDeform_Mesh()) {
SU2_MPI::Error("Initial coordinates are only available with DEFORM_MESH= YES", CURRENT_FUNCTION);
}
auto* coords =
const_cast<su2activematrix*>(solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord());
return CPyWrapperMatrixView(*coords, "InitialCoordinates", true);
}

/*!
* \brief Set the coordinates of a mesh node.
* \param[in] iPoint - Mesh node index.
* \param[in] values - Node coordinates (nDim).
* \brief Get a read/write view of the current coordinates of all mesh nodes.
*/
void SetCoordinates(unsigned long iPoint, vector<passivedouble> values);
inline CPyWrapperMatrixView Coordinates() {
auto& coords = const_cast<su2activematrix&>(main_geometry->nodes->GetCoord());
return CPyWrapperMatrixView(coords, "Coordinates", false);
}

/*!
* \brief Get the number of markers in the mesh.
Expand Down Expand Up @@ -323,6 +322,39 @@ class CDriverBase {
* \brief Communicate the boundary mesh displacements.
*/
void CommunicateMeshDisplacements(void);

/*!
* \brief Get all the active solver names with their associated indices (which can be used to access their data).
*/
map<string, unsigned short> GetSolverIndices() const;

/*!
* \brief Get a read/write view of the current solution on all mesh nodes of a solver.
*/
inline CPyWrapperMatrixView Solution(unsigned short iSolver) {
auto* solver = solver_container[ZONE_0][INST_0][MESH_0][iSolver];
if (solver == nullptr) SU2_MPI::Error("The selected solver does not exist.", CURRENT_FUNCTION);
auto& solution = solver->GetNodes()->GetSolution();
return CPyWrapperMatrixView(solution, "Solution of " + solver->GetSolverName(), false);
}

/*!
* \brief Get the flow solver primitive variable names with their associated indices.
* These correspond to the column indices in the matrix returned by Primitives.
*/
map<string, unsigned short> GetPrimitiveIndices() const;

/*!
* \brief Get a read/write view of the current primitive variables on all mesh nodes of the flow solver.
* \warning Primitive variables are only available for flow solvers.
*/
inline CPyWrapperMatrixView Primitives() {
auto* solver = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL];
if (solver == nullptr) SU2_MPI::Error("The flow solver does not exist.", CURRENT_FUNCTION);
auto& primitives = const_cast<su2activematrix&>(solver->GetNodes()->GetPrimitive());
return CPyWrapperMatrixView(primitives, "Primitives", false);
}

/// \}

protected:
Expand Down
70 changes: 5 additions & 65 deletions SU2_CFD/include/output/CFlowOutput.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,10 @@
#include "CFVMOutput.hpp"
#include "../variables/CVariable.hpp"

/*--- Forward declare to avoid including here. ---*/
template <class>
struct CPrimitiveIndices;

class CFlowOutput : public CFVMOutput{
protected:
unsigned long lastInnerIter;
Expand Down Expand Up @@ -211,71 +215,7 @@ class CFlowOutput : public CFVMOutput{
* \brief Helper for custom outputs, converts variable names to indices and pointers which are then used
* to evaluate the custom expressions.
*/
template <class FlowIndices>
void ConvertVariableSymbolsToIndices(const FlowIndices& idx, CustomOutput& output) const {

static const auto knownVariables =
"TEMPERATURE, TEMPERATURE_VE, VELOCITY_X, VELOCITY_Y, VELOCITY_Z, PRESSURE,\n"
"DENSITY, ENTHALPY, SOUND_SPEED, LAMINAR_VISCOSITY, EDDY_VISCOSITY, THERMAL_CONDUCTIVITY\n"
"TURB[0,1,...], RAD[0,1,...], SPECIES[0,1,...]";

auto IndexOfVariable = [&](const FlowIndices& idx, const std::string& var) {
/*--- Primitives of the flow solver. ---*/
const auto flow_offset = FLOW_SOL * CustomOutput::MAX_VARS_PER_SOLVER;

if ("TEMPERATURE" == var) return flow_offset + idx.Temperature();
if ("TEMPERATURE_VE" == var) return flow_offset + idx.Temperature_ve();
if ("VELOCITY_X" == var) return flow_offset + idx.Velocity();
if ("VELOCITY_Y" == var) return flow_offset + idx.Velocity() + 1;
if ("VELOCITY_Z" == var) return flow_offset + idx.Velocity() + 2;
if ("PRESSURE" == var) return flow_offset + idx.Pressure();
if ("DENSITY" == var) return flow_offset + idx.Density();
if ("ENTHALPY" == var) return flow_offset + idx.Enthalpy();
if ("SOUND_SPEED" == var) return flow_offset + idx.SoundSpeed();
if ("LAMINAR_VISCOSITY" == var) return flow_offset + idx.LaminarViscosity();
if ("EDDY_VISCOSITY" == var) return flow_offset + idx.EddyViscosity();
if ("THERMAL_CONDUCTIVITY" == var) return flow_offset + idx.ThermalConductivity();

/*--- Index-based (no name) access to variables of other solvers. ---*/
auto GetIndex = [](const std::string& s, int nameLen) {
/*--- Extract an int from "name[int]", nameLen is the length of "name". ---*/
return std::stoi(std::string(s.begin() + nameLen + 1, s.end() - 1));
};
if (var.rfind("SPECIES", 0) == 0) return SPECIES_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 7);
if (var.rfind("TURB", 0) == 0) return TURB_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 4);
if (var.rfind("RAD", 0) == 0) return RAD_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 3);

return CustomOutput::NOT_A_VARIABLE;
};

output.otherOutputs.clear();
output.varIndices.clear();
output.varIndices.reserve(output.varSymbols.size());

for (const auto& var : output.varSymbols) {
output.varIndices.push_back(IndexOfVariable(idx, var));

if (output.type == OperationType::FUNCTION && output.varIndices.back() != CustomOutput::NOT_A_VARIABLE) {
SU2_MPI::Error("Custom outputs of type 'Function' cannot reference solver variables.", CURRENT_FUNCTION);
}
/*--- Symbol is a valid solver variable. ---*/
if (output.varIndices.back() < CustomOutput::NOT_A_VARIABLE) continue;

/*--- An index above NOT_A_VARIABLE is not valid with current solver settings. ---*/
if (output.varIndices.back() > CustomOutput::NOT_A_VARIABLE) {
SU2_MPI::Error("Inactive solver variable (" + var + ") used in function " + output.name + "\n"
"E.g. this may only be a variable of the compressible solver.", CURRENT_FUNCTION);
}

/*--- An index equal to NOT_A_VARIABLE may refer to a history output. ---*/
output.varIndices.back() += output.otherOutputs.size();
output.otherOutputs.push_back(GetPtrToHistoryOutput(var));
if (output.otherOutputs.back() == nullptr) {
SU2_MPI::Error("Invalid history output or solver variable (" + var + ") used in function " + output.name +
"\nValid solvers variables: " + knownVariables, CURRENT_FUNCTION);
}
}
}
void ConvertVariableSymbolsToIndices(const CPrimitiveIndices<unsigned long>& idx, CustomOutput& output) const;

/*!
* \brief Compute value of the Q criteration for vortex idenfitication
Expand Down
4 changes: 2 additions & 2 deletions SU2_CFD/include/solvers/CSolver.hpp
View file
Open in desktop
Loading