[WIP] Modernize scalar numerics

Common/include/containers/container_decorators.hpp

@@ -62,6 +62,33 @@ class CMatrixView {
   const Scalar* operator[](Index i) const noexcept { return &m_ptr[i * m_cols]; }
   const Scalar& operator()(Index i, Index j) const noexcept { return m_ptr[i * m_cols + j]; }
 
+  /*!
+   * \brief Return copy of data in a static size container (see C2DContainer::get).
+   * \param[in] i - Row of the view (e.g. point index, whole-mesh usage).
+   * \param[in] start - Starting column to copy the data (amount determined by container size).
+   */
+  template <class StaticContainer>
+  StaticContainer get(Index i, Index start = 0) const noexcept {
+    constexpr size_t Size = StaticContainer::StaticSize;
+    static_assert(Size, "This method requires a static output type.");
+    StaticContainer ret;
+    for (size_t k = 0; k < Size; ++k) ret.data()[k] = m_ptr[i * m_cols + start + k];
+    return ret;
+  }
+
+  /*!
+   * \brief SIMD gather version of get, one row per lane.
+   */
+  template <class StaticContainer, class U, size_t N>
+  StaticContainer get(simd::Array<U, N> i, Index start = 0) const noexcept {
+    constexpr size_t Size = StaticContainer::StaticSize;
+    static_assert(Size, "This method requires a static output type.");
+    StaticContainer ret;
+    for (size_t lane = 0; lane < N; ++lane)
+      for (size_t k = 0; k < Size; ++k) ret.data()[k][lane] = m_ptr[i[lane] * m_cols + start + k];
+    return ret;
+  }
+
   template <class U = T, su2enable_if<!std::is_const<U>::value> = 0>
   Scalar* operator[](Index i) noexcept {
     return &m_ptr[i * m_cols];

Common/include/linear_algebra/CSysMatrix.hpp

@@ -1027,6 +1027,80 @@ class CSysMatrix {
     SetBlocks<MatrixType, ScalarType, false>(iEdge, block_i, block_j, -1);
   }
 
+  /*!
+   * \brief Set the four blocks of an edge, for fluxes whose i and j contributions are independent.
+   * \note The diagonal blocks are accumulated, the off-diagonal blocks are set.
+   */
+  template <class MatrixType, class OtherType = ScalarType>
+  inline void SetBlocks(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, const MatrixType& jac_ii,
+                        const MatrixType& jac_ij, const MatrixType& jac_ji, const MatrixType& jac_jj,
+                        OtherType mask = 1) {
+    const auto blkSz = nVar * nEqn;
+    auto* bii = &mat.d[iPoint * blkSz];
+    auto* bjj = &mat.d[jPoint * blkSz];
+    unsigned long iVar, jVar, offset = 0;
+
+    if (quantized_mode) {
+      ScalarType bij_buf[MAXNVAR * MAXNVAR], bji_buf[MAXNVAR * MAXNVAR];
+      for (iVar = 0; iVar < nVar; iVar++)
+        for (jVar = 0; jVar < nEqn; jVar++, ++offset) {
+          bii[offset] += PassiveAssign(jac_ii[iVar][jVar] * mask);
+          bjj[offset] += PassiveAssign(jac_jj[iVar][jVar] * mask);
+          bij_buf[offset] = PassiveAssign(jac_ij[iVar][jVar] * mask);
+          bji_buf[offset] = PassiveAssign(jac_ji[iVar][jVar] * mask);
+        }
+      QuantizeBlock(bij_buf, &q_scale.u[iEdge * nVar], &q_blocks.u[iEdge * blkSz]);
+      const auto k_l = edge_ptr_l[iEdge];
+      QuantizeBlock(bji_buf, &q_scale.l[k_l * nVar], &q_blocks.l[k_l * blkSz]);
+      return;
+    }
+
+    auto* bij = &mat.u[iEdge * blkSz];
+    auto* bji = &mat.l[edge_ptr_l[iEdge] * blkSz];
+    for (iVar = 0; iVar < nVar; iVar++) {
+      for (jVar = 0; jVar < nEqn; jVar++) {
+        bii[offset] += PassiveAssign(jac_ii[iVar][jVar] * mask);
+        bjj[offset] += PassiveAssign(jac_jj[iVar][jVar] * mask);
+        bij[offset] = PassiveAssign(jac_ij[iVar][jVar] * mask);
+        bji[offset] = PassiveAssign(jac_ji[iVar][jVar] * mask);
+        ++offset;
+      }
+    }
+  }
+
+  /*!
+   * \brief Set the off-diagonal blocks of an edge, the diagonal being assembled elsewhere.
+   */
+  template <class MatrixType, class OtherType = ScalarType>
+  inline void SetOffDiagBlocks(unsigned long iEdge, const MatrixType& jac_ij, const MatrixType& jac_ji,
+                               OtherType mask = 1) {
+    const auto blkSz = nVar * nEqn;
+    unsigned long iVar, jVar, offset = 0;
+
+    if (quantized_mode) {
+      ScalarType bij_buf[MAXNVAR * MAXNVAR], bji_buf[MAXNVAR * MAXNVAR];
+      for (iVar = 0; iVar < nVar; iVar++)
+        for (jVar = 0; jVar < nEqn; jVar++, ++offset) {
+          bij_buf[offset] = PassiveAssign(jac_ij[iVar][jVar] * mask);
+          bji_buf[offset] = PassiveAssign(jac_ji[iVar][jVar] * mask);
+        }
+      QuantizeBlock(bij_buf, &q_scale.u[iEdge * nVar], &q_blocks.u[iEdge * blkSz]);
+      const auto k_l = edge_ptr_l[iEdge];
+      QuantizeBlock(bji_buf, &q_scale.l[k_l * nVar], &q_blocks.l[k_l * blkSz]);
+      return;
+    }
+
+    auto* bij = &mat.u[iEdge * blkSz];
+    auto* bji = &mat.l[edge_ptr_l[iEdge] * blkSz];
+    for (iVar = 0; iVar < nVar; iVar++) {
+      for (jVar = 0; jVar < nEqn; jVar++) {
+        bij[offset] = PassiveAssign(jac_ij[iVar][jVar] * mask);
+        bji[offset] = PassiveAssign(jac_ji[iVar][jVar] * mask);
+        ++offset;
+      }
+    }
+  }
+
   /*!
    * \brief SIMD version, does the update for multiple edges.
    * \note Nothing is updated if the mask is 0.

Common/src/CConfig.cpp

@@ -4352,6 +4352,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     }
 
     /* --- Check for NEMO compatibility issues ---*/
+    if (nemo && Kind_Turb_Model != TURB_MODEL::NONE) {
+      SU2_MPI::Error("A turbulence model is not yet available for the NEMO solver.", CURRENT_FUNCTION);
+    }
+
     if (Kind_FluidModel == SU2_NONEQ && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::SUTHERLAND && Kind_TransCoeffModel != TRANSCOEFFMODEL::GUPTAYOS) ) {
       SU2_MPI::Error("Transport model not available for NEMO solver using SU2TCLIB. Please use the WILKE, SUTHERLAND or GUPTAYOS transport model instead.", CURRENT_FUNCTION);
     }