Newton-Krylov primal iterations & Krylov discrete adjoint

Common/include/CConfig.hpp

@@ -28,7 +28,7 @@
 
 #pragma once
 
-#include "./parallelization/mpi_structure.hpp"
+#include "parallelization/mpi_structure.hpp"
 
 #include <iostream>
 #include <cstdlib>
@@ -416,6 +416,9 @@ class CConfig {
 
   unsigned short nQuasiNewtonSamples;  /*!< \brief Number of samples used in quasi-Newton solution methods. */
   bool UseVectorization;       /*!< \brief Whether to use vectorized numerics schemes. */
+  bool NewtonKrylov;           /*!< \brief Use a coupled Newton method to solve the flow equations. */
+  array<unsigned short,3> NK_IntParam{{20, 3, 2}}; /*!< \brief Integer parameters for NK method. */
+  array<su2double,4> NK_DblParam{{-2.0, 0.1, -3.0, 1e-4}}; /*!< \brief Floating-point parameters for NK method. */
 
   unsigned short nMGLevels;    /*!< \brief Number of multigrid levels (coarse levels). */
   unsigned short nCFL;         /*!< \brief Number of CFL, one for each multigrid level. */
@@ -1185,6 +1188,8 @@ class CConfig {
 
   void addDoubleArrayOption(const string name, const int size, su2double* option_field);
 
+  void addUShortArrayOption(const string name, const int size, unsigned short* option_field);
+
   void addDoubleListOption(const string name, unsigned short & size, su2double * & option_field);
 
   void addShortListOption(const string name, unsigned short & size, short * & option_field);
@@ -3970,6 +3975,21 @@ class CConfig {
    */
   bool GetUseVectorization(void) const { return UseVectorization; }
 
+  /*!
+   * \brief Get whether to use a Newton-Krylov method.
+   */
+  bool GetNewtonKrylov(void) const { return NewtonKrylov; }
+
+  /*!
+   * \brief Get Newton-Krylov integer parameters.
+   */
+  array<unsigned short,3> GetNewtonKrylovIntParam(void) const { return NK_IntParam; }
+
+  /*!
+   * \brief Get Newton-Krylov floating-point parameters.
+   */
+  array<su2double,4> GetNewtonKrylovDblParam(void) const { return NK_DblParam; }
+
   /*!
    * \brief Get the relaxation coefficient of the linear solver for the implicit formulation.
    * \return relaxation coefficient of the linear solver for the implicit formulation.

Common/include/linear_algebra/CSysMatrix.hpp

@@ -28,9 +28,7 @@
 
 #pragma once
 
-#include "../../include/parallelization/mpi_structure.hpp"
-#include "../../include/parallelization/omp_structure.hpp"
-#include "../../include/parallelization/vectorization.hpp"
+#include "../../include/CConfig.hpp"
 #include "CSysVector.hpp"
 #include "CPastixWrapper.hpp"
 
@@ -75,16 +73,43 @@ struct mkl_jit_wrapper<float> {
 #endif
 #endif
 
-class CConfig;
 class CGeometry;
 
+struct CSysMatrixComms {
+  /*!
+   * \brief Routine to load a vector quantity into the data structures for MPI point-to-point
+   *        communication and to launch non-blocking sends and recvs.
+   * \param[in] x        - CSysVector holding the array of data.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config   - Definition of the particular problem.
+   * \param[in] commType - Enumerated type for the quantity to be communicated.
+   */
+  template<class T>
+  static void Initiate(const CSysVector<T>& x, CGeometry *geometry, const CConfig *config,
+                       unsigned short commType = SOLUTION_MATRIX);
+
+  /*!
+   * \brief Routine to complete the set of non-blocking communications launched by
+   *        Initiate() and unpacking of the data in the vector.
+   * \param[in] x        - CSysVector holding the array of data.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config   - Definition of the particular problem.
+   * \param[in] commType - Enumerated type for the quantity to be unpacked.
+   */
+  template<class T>
+  static void Complete(CSysVector<T>& x, CGeometry *geometry, const CConfig *config,
+                       unsigned short commType = SOLUTION_MATRIX);
+};
+
 /*!
  * \class CSysMatrix
  * \brief Main class for defining block-compressed-row-storage sparse matrices.
  */
 template<class ScalarType>
 class CSysMatrix {
 private:
+  friend class CSysMatrixComms;
+
   const int rank;     /*!< \brief MPI Rank. */
   const int size;     /*!< \brief MPI Size. */
 
@@ -163,9 +188,11 @@ class CSysMatrix {
 
   /*!
    * \brief Handle type conversion for when we Set, Add, etc. blocks, preserving derivative information (if supported by types).
-   * \note See specialization for discrete adjoint right outside this class's declaration.
    */
-  template<class DstType, class SrcType>
+  template<class DstType, class SrcType, su2enable_if<std::is_arithmetic<DstType>::value> = 0>
+  FORCEINLINE static DstType ActiveAssign(const SrcType& val) { return SU2_TYPE::GetValue(val); }
+
+  template<class DstType, class SrcType, su2enable_if<!std::is_arithmetic<DstType>::value> = 0>
   FORCEINLINE static DstType ActiveAssign(const SrcType& val) { return val; }
 
   /*!
@@ -378,34 +405,6 @@ class CSysMatrix {
    */
   void SetValDiagonalZero(void);
 
-  /*!
-   * \brief Routine to load a vector quantity into the data structures for MPI point-to-point
-   *        communication and to launch non-blocking sends and recvs.
-   * \param[in] x        - CSysVector holding the array of data.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config   - Definition of the particular problem.
-   * \param[in] commType - Enumerated type for the quantity to be communicated.
-   */
-  template<class OtherType>
-  void InitiateComms(const CSysVector<OtherType> & x,
-                     CGeometry *geometry,
-                     const CConfig *config,
-                     unsigned short commType) const;
-
-  /*!
-   * \brief Routine to complete the set of non-blocking communications launched by
-   *        InitiateComms() and unpacking of the data in the vector.
-   * \param[in] x        - CSysVector holding the array of data.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config   - Definition of the particular problem.
-   * \param[in] commType - Enumerated type for the quantity to be unpacked.
-   */
-  template<class OtherType>
-  void CompleteComms(CSysVector<OtherType> & x,
-                     CGeometry *geometry,
-                     const CConfig *config,
-                     unsigned short commType) const;
-
   /*!
    * \brief Get a pointer to the start of block "ij"
    * \param[in] block_i - Row index.
@@ -918,8 +917,3 @@ class CSysMatrix {
                                    CGeometry *geometry, const CConfig *config) const;
 
 };
-
-#ifdef CODI_REVERSE_TYPE
-template<> template<>
-FORCEINLINE su2mixedfloat CSysMatrix<su2mixedfloat>::ActiveAssign(const su2double& val) { return SU2_TYPE::GetValue(val); }
-#endif

Common/include/linear_algebra/CSysSolve.hpp

@@ -81,7 +81,6 @@ class CSysSolve {
 
   mutable bool cg_ready;     /*!< \brief Indicate if memory used by CG is allocated. */
   mutable bool bcg_ready;    /*!< \brief Indicate if memory used by BCGSTAB is allocated. */
-  mutable bool gmres_ready;  /*!< \brief Indicate if memory used by FGMRES is allocated. */
   mutable bool smooth_ready; /*!< \brief Indicate if memory used by SMOOTHER is allocated. */
 
   mutable VectorType r;      /*!< \brief Residual in CG and BCGSTAB. */
@@ -101,6 +100,9 @@ class CSysSolve {
   const VectorType* LinSysRes_ptr;  /*!< \brief Pointer to appropriate LinSysRes (set to original or temporary in call to Solve). */
 
   LinearToleranceType tol_type = LinearToleranceType::RELATIVE; /*!< \brief How the linear solvers interpret the tolerance. */
+  bool xIsZero = false;           /*!< \brief If true assume the initial solution is always 0. */
+  bool recomputeRes = false;      /*!< \brief Recompute the residual after inner iterations, if monitoring. */
+  unsigned long monitorFreq = 10; /*!< \brief Monitoring frequency. */
 
   /*!
    * \brief sign transfer function
@@ -388,4 +390,19 @@ class CSysSolve {
    */
   inline void SetToleranceType(LinearToleranceType type) {tol_type = type;}
 
+  /*!
+   * \brief Assume the initial solution is 0 to save one product, or don't.
+   */
+  inline void SetxIsZero(bool isZero) {xIsZero = isZero;}
+
+  /*!
+   * \brief Set whether to recompute residuals at the end (while monitoring only).
+   */
+  inline void SetRecomputeResidual(bool recompRes) {recomputeRes = recompRes;}
+
+  /*!
+   * \brief Set the screen output frequency during monitoring.
+   */
+  inline void SetMonitoringFrequency(bool frequency) {monitorFreq = frequency;}
+
 };

Common/include/linear_algebra/CSysVector.hpp

@@ -66,7 +66,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
   ScalarType* vec_val = nullptr;               /*!< \brief Storage, 64 byte aligned (do not use normal new/delete). */
   unsigned long nElm = 0;          /*!< \brief Total number of elements (or number elements on this processor). */
   unsigned long nElmDomain = 0;    /*!< \brief Total number of elements without Ghost cells. */
-  unsigned long nVar = 0;          /*!< \brief Number of elements in a block. */
+  unsigned long nVar = 1;          /*!< \brief Number of elements in a block. */
 
   /*!
    * \brief Generic initialization from a scalar or array.
@@ -111,7 +111,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \param[in] size - Number of elements locally.
    * \param[in] val - Default value for elements.
    */
-  CSysVector(unsigned long size, ScalarType val = 0.0) { Initialize(size, size, 1, &val, false); }
+  explicit CSysVector(unsigned long size, ScalarType val = 0.0) { Initialize(size, size, 1, &val, false); }
 
   /*!
    * \brief Construct from size and value (block version).
@@ -129,7 +129,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \param[in] size - Number of elements locally.
    * \param[in] u_array - Vector stored as array being copied.
    */
-  explicit CSysVector(unsigned long size, const ScalarType* u_array) { Initialize(size, size, 1, u_array, true); }
+  CSysVector(unsigned long size, const ScalarType* u_array) { Initialize(size, size, 1, u_array, true); }
 
   /*!
    * \brief Constructor from array (block version).
@@ -138,8 +138,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \param[in] numVar - number of variables in each block
    * \param[in] u_array - vector stored as array being copied
    */
-  explicit CSysVector(unsigned long numBlk, unsigned long numBlkDomain, unsigned long numVar,
-                      const ScalarType* u_array) {
+  CSysVector(unsigned long numBlk, unsigned long numBlkDomain, unsigned long numVar, const ScalarType* u_array) {
     Initialize(numBlk, numBlkDomain, numVar, u_array, true);
   }
 

Common/include/option_structure.hpp

@@ -77,7 +77,7 @@ const unsigned int MAX_PARAMETERS = 10;       /*!< \brief Maximum number of para
 const unsigned int MAX_NUMBER_PERIODIC = 10;  /*!< \brief Maximum number of periodic boundary conditions. */
 const unsigned int MAX_STRING_SIZE = 200;     /*!< \brief Maximum number of domains. */
 const unsigned int MAX_NUMBER_FFD = 15;       /*!< \brief Maximum number of FFDBoxes for the FFD. */
-const unsigned int MAX_SOLS = 12;             /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */
+enum: unsigned int{MAX_SOLS = 12};            /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */
 const unsigned int MAX_TERMS = 6;             /*!< \brief Maximum number of terms in the numerical equations (dimension of solver container array). */
 const unsigned int MAX_ZONES = 3;             /*!< \brief Maximum number of zones. */
 const unsigned int MAX_FE_KINDS = 4;          /*!< \brief Maximum number of Finite Elements. */
@@ -235,7 +235,7 @@ static const MapType<string, ENUM_MAIN_SOLVER> Solver_Map = {
 };
 
 /*!
- * \brief different solver types for the multizone environment component
+ * \brief Different solver types for multizone problems
  */
 enum ENUM_MULTIZONE {
   MZ_BLOCK_GAUSS_SEIDEL = 0,   /*!< \brief Definition of a Block-Gauss-Seidel multizone solver. */
@@ -439,7 +439,6 @@ static const MapType<string, ENUM_MEASUREMENTS> Measurements_Map = {
 enum RUNTIME_TYPE {
   RUNTIME_FLOW_SYS = 2,       /*!< \brief One-physics case, the code is solving the flow equations(Euler and Navier-Stokes). */
   RUNTIME_TURB_SYS = 3,       /*!< \brief One-physics case, the code is solving the turbulence model. */
-  RUNTIME_ADJPOT_SYS = 5,     /*!< \brief One-physics case, the code is solving the adjoint potential flow equation. */
   RUNTIME_ADJFLOW_SYS = 6,    /*!< \brief One-physics case, the code is solving the adjoint equations is being solved (Euler and Navier-Stokes). */
   RUNTIME_ADJTURB_SYS = 7,    /*!< \brief One-physics case, the code is solving the adjoint turbulence model. */
   RUNTIME_MULTIGRID_SYS = 14, /*!< \brief Full Approximation Storage Multigrid system of equations. */

Common/include/option_structure.inl

@@ -334,19 +334,20 @@ public:
   }
 };
 
-class COptionDoubleArray : public COptionBase {
+template<class Type>
+class COptionArray : public COptionBase {
   string name; // Identifier for the option
   const int size; // Number of elements
-  su2double* field; // Reference to the fieldname
+  Type* field; // Reference to the field
 
 public:
-  COptionDoubleArray(string option_field_name, const int list_size, su2double* option_field) :
+  COptionArray(string option_field_name, const int list_size, Type* option_field) :
     name(option_field_name),
     size(list_size),
     field(option_field) {
   }
 
-  ~COptionDoubleArray() override {};
+  ~COptionArray() override {};
 
   string SetValue(vector<string> option_value) override {
     COptionBase::SetValue(option_value);
@@ -368,7 +369,7 @@ public:
     for (int i  = 0; i < this->size; i++) {
       istringstream is(option_value[i]);
       if (!(is >> field[i])) {
-        return badValue(option_value, "su2double array", this->name);
+        return badValue(option_value, " array", this->name);
       }
     }
     return "";

Common/include/parallelization/mpi_structure.cpp

@@ -122,6 +122,49 @@ void CBaseMPIWrapper::Error(std::string ErrorMsg, std::string FunctionName){
   }
   Abort(currentComm, 0);
 }
+
+void CBaseMPIWrapper::CopyData(const void* sendbuf, void* recvbuf, int size, Datatype datatype) {
+  switch (datatype) {
+    case MPI_DOUBLE:
+      for (int i = 0; i < size; i++) {
+        static_cast<su2double*>(recvbuf)[i] = static_cast<const su2double*>(sendbuf)[i];
+      }
+      break;
+    case MPI_UNSIGNED_LONG:
+      for (int i = 0; i < size; i++) {
+        static_cast<unsigned long*>(recvbuf)[i] = static_cast<const unsigned long*>(sendbuf)[i];
+      }
+      break;
+    case MPI_LONG:
+      for (int i = 0; i < size; i++) {
+        static_cast<long*>(recvbuf)[i] = static_cast<const long*>(sendbuf)[i];
+      }
+      break;
+    case MPI_UNSIGNED_SHORT:
+      for (int i = 0; i < size; i++) {
+        static_cast<unsigned short*>(recvbuf)[i] = static_cast<const unsigned short*>(sendbuf)[i];
+      }
+      break;
+    case MPI_CHAR:
+      for (int i = 0; i < size; i++) {
+        static_cast<char*>(recvbuf)[i] = static_cast<const char*>(sendbuf)[i];
+      }
+      break;
+    case MPI_SHORT:
+      for (int i = 0; i < size; i++) {
+        static_cast<short*>(recvbuf)[i] = static_cast<const short*>(sendbuf)[i];
+      }
+      break;
+    case MPI_INT:
+      for (int i = 0; i < size; i++) {
+        static_cast<int*>(recvbuf)[i] = static_cast<const int*>(sendbuf)[i];
+      }
+      break;
+    default:
+      Error("Unknown type", CURRENT_FUNCTION);
+      break;
+  };
+}
 #endif
 
 #ifdef HAVE_MPI

Common/include/parallelization/mpi_structure.hpp

@@ -503,48 +503,7 @@ class CBaseMPIWrapper {
   static int Rank, Size;
   static Comm currentComm;
 
-  static inline void CopyData(const void* sendbuf, void* recvbuf, int size, Datatype datatype) {
-    switch (datatype) {
-      case MPI_DOUBLE:
-        for (int i = 0; i < size; i++) {
-          static_cast<su2double*>(recvbuf)[i] = static_cast<const su2double*>(sendbuf)[i];
-        }
-        break;
-      case MPI_UNSIGNED_LONG:
-        for (int i = 0; i < size; i++) {
-          static_cast<unsigned long*>(recvbuf)[i] = static_cast<const unsigned long*>(sendbuf)[i];
-        }
-        break;
-      case MPI_LONG:
-        for (int i = 0; i < size; i++) {
-          static_cast<long*>(recvbuf)[i] = static_cast<const long*>(sendbuf)[i];
-        }
-        break;
-      case MPI_UNSIGNED_SHORT:
-        for (int i = 0; i < size; i++) {
-          static_cast<unsigned short*>(recvbuf)[i] = static_cast<const unsigned short*>(sendbuf)[i];
-        }
-        break;
-      case MPI_CHAR:
-        for (int i = 0; i < size; i++) {
-          static_cast<char*>(recvbuf)[i] = static_cast<const char*>(sendbuf)[i];
-        }
-        break;
-      case MPI_SHORT:
-        for (int i = 0; i < size; i++) {
-          static_cast<short*>(recvbuf)[i] = static_cast<const short*>(sendbuf)[i];
-        }
-        break;
-      case MPI_INT:
-        for (int i = 0; i < size; i++) {
-          static_cast<int*>(recvbuf)[i] = static_cast<const int*>(sendbuf)[i];
-        }
-        break;
-      default:
-        Error("Unknown type", CURRENT_FUNCTION);
-        break;
-    };
-  }
+  static void CopyData(const void* sendbuf, void* recvbuf, int size, Datatype datatype);
 
  public:
   static void Error(std::string ErrorMsg, std::string FunctionName);