Access solution and primitives via python wrapper (#1938)

* matrix view type to manipulate fields in python without copying data

* extra line and update docker images

Author: pcarruscag

.github/workflows/regression.yml

@@ -52,12 +52,12 @@ jobs:
           key: ${{ matrix.config_set }}-${{ github.sha }}
           restore-keys: ${{ matrix.config_set }}
       - name: Pre Cleanup
-        uses: docker://ghcr.io/su2code/su2/build-su2:221224-1158
+        uses: docker://ghcr.io/su2code/su2/build-su2:230225-2136
         with:
           entrypoint: /bin/rm
           args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}
       - name: Build
-        uses: docker://ghcr.io/su2code/su2/build-su2:221224-1158
+        uses: docker://ghcr.io/su2code/su2/build-su2:230225-2136
         with:
           args: -b ${{github.ref}} -f "${{matrix.flags}}"
       - name: Compress binaries
@@ -68,7 +68,7 @@ jobs:
           name: ${{ matrix.config_set }}
           path: install_bin.tgz
       - name: Post Cleanup
-        uses: docker://ghcr.io/su2code/su2/build-su2:221224-1158
+        uses: docker://ghcr.io/su2code/su2/build-su2:230225-2136
         with:
           entrypoint: /bin/rm
           args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}
@@ -99,7 +99,7 @@ jobs:
             tag: OMP
     steps:
       - name: Pre Cleanup
-        uses: docker://ghcr.io/su2code/su2/test-su2:221224-1158
+        uses: docker://ghcr.io/su2code/su2/test-su2:230225-2136
         with:
           entrypoint: /bin/rm
           args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}
@@ -125,12 +125,12 @@ jobs:
           chmod a+x $BIN_FOLDER/*
           ls -lahR $BIN_FOLDER
       - name: Run Tests in Container
-        uses: docker://ghcr.io/su2code/su2/test-su2:221224-1158
+        uses: docker://ghcr.io/su2code/su2/test-su2:230225-2136
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
           args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
       - name: Cleanup
-        uses: docker://ghcr.io/su2code/su2/test-su2:221224-1158
+        uses: docker://ghcr.io/su2code/su2/test-su2:230225-2136
         with:
           entrypoint: /bin/rm
           args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}
@@ -151,7 +151,7 @@ jobs:
             tag: MPI
     steps:
       - name: Pre Cleanup
-        uses: docker://ghcr.io/su2code/su2/test-su2:221224-1158
+        uses: docker://ghcr.io/su2code/su2/test-su2:230225-2136
         with:
           entrypoint: /bin/rm
           args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}
@@ -212,11 +212,11 @@ jobs:
           echo $PWD
           ls -lahR
       - name: Run Unit Tests
-        uses: docker://ghcr.io/su2code/su2/test-su2:221224-1158
+        uses: docker://ghcr.io/su2code/su2/test-su2:230225-2136
         with:
           entrypoint: install/bin/${{matrix.testdriver}}
       - name: Post Cleanup
-        uses: docker://ghcr.io/su2code/su2/test-su2:221224-1158
+        uses: docker://ghcr.io/su2code/su2/test-su2:230225-2136
         with:
           entrypoint: /bin/rm
           args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}

.github/workflows/release-management.yml

@@ -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

Common/include/containers/CPyWrapperMatrixView.hpp

@@ -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);
+  }
+};

Common/src/containers/meson.build

@@ -1,3 +1,3 @@
 common_src += files(['CTrapezoidalMap.cpp',
-                     'CFileReaderLUT.cpp', 
+                     'CFileReaderLUT.cpp',
                      'CLookUpTable.cpp'])

SU2_CFD/include/drivers/CDriverBase.hpp

@@ -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"
@@ -38,7 +39,6 @@
  * \brief Base class for all drivers.
  * \author H. Patel, A. Gastaldi
  */
-
 class CDriverBase {
  protected:
   int rank,               /*!< \brief MPI Rank. */
@@ -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.
@@ -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:

SU2_CFD/include/output/CFlowOutput.hpp

@@ -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;
@@ -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

SU2_CFD/include/solvers/CSolver.hpp

@@ -4246,9 +4246,9 @@ class CSolver {
 
   /*!
    * \brief Retrieve the solver name for output purposes.
-   * \param[out] val_solvername - Name of the solver.
+   * \returns Name of the solver.
    */
-  inline string GetSolverName(void) {return SolverName;}
+  inline const string& GetSolverName() const  { return SolverName; }
 
   /*!
    * \brief Get the solution fields.