diff --git a/Common/include/config_structure.hpp b/Common/include/config_structure.hpp index e2703516568c..586a856923f9 100644 --- a/Common/include/config_structure.hpp +++ b/Common/include/config_structure.hpp @@ -554,6 +554,7 @@ class CConfig { su2double Deform_Linear_Solver_Error; /*!< \brief Min error of the linear solver for the implicit formulation. */ su2double Linear_Solver_Error_FSI_Struc; /*!< \brief Min error of the linear solver for the implicit formulation in the structural side for FSI problems . */ su2double Linear_Solver_Error_Heat; /*!< \brief Min error of the linear solver for the implicit formulation in the fvm heat solver . */ + su2double Linear_Solver_Smoother_Relaxation; /*!< \brief Relaxation factor for iterative linear smoothers. */ unsigned long Linear_Solver_Iter; /*!< \brief Max iterations of the linear solver for the implicit formulation. */ unsigned long Deform_Linear_Solver_Iter; /*!< \brief Max iterations of the linear solver for the implicit formulation. */ unsigned long Linear_Solver_Iter_FSI_Struc; /*!< \brief Max iterations of the linear solver for FSI applications and structural solver. */ @@ -4025,6 +4026,12 @@ class CConfig { */ unsigned long GetLinear_Solver_Restart_Frequency(void); + /*! + * \brief Get the relaxation factor for iterative linear smoothers. + * \return Relaxation factor. + */ + su2double GetLinear_Solver_Smoother_Relaxation(void) const; + /*! * \brief Get the relaxation coefficient of the linear solver for the implicit formulation. * \return relaxation coefficient of the linear solver for the implicit formulation. diff --git a/Common/include/config_structure.inl b/Common/include/config_structure.inl index 5ed8e1f2c0e2..08c21294fe9f 100644 --- a/Common/include/config_structure.inl +++ b/Common/include/config_structure.inl @@ -1035,6 +1035,8 @@ inline unsigned short CConfig::GetLinear_Solver_ILU_n(void) { return Linear_Solv inline unsigned long CConfig::GetLinear_Solver_Restart_Frequency(void) { return Linear_Solver_Restart_Frequency; } +inline su2double CConfig::GetLinear_Solver_Smoother_Relaxation(void) const { return Linear_Solver_Smoother_Relaxation; } + inline su2double CConfig::GetRelaxation_Factor_Flow(void) { return Relaxation_Factor_Flow; } inline su2double CConfig::GetRelaxation_Factor_AdjFlow(void) { return Relaxation_Factor_AdjFlow; } diff --git a/Common/include/interpolation_structure.hpp b/Common/include/interpolation_structure.hpp index 0e937bd93422..7149262701b9 100644 --- a/Common/include/interpolation_structure.hpp +++ b/Common/include/interpolation_structure.hpp @@ -50,14 +50,6 @@ #include "geometry_structure.hpp" #include "vector_structure.hpp" -#ifdef HAVE_LAPACK -/*--- Lapack / Blas routines used in RBF interpolation. ---*/ -extern "C" void dsptrf_(char*, int*, passivedouble*, int*, int*); -extern "C" void dsptri_(char*, int*, passivedouble*, int*, passivedouble*, int*); -extern "C" void dsymm_(char*, char*, int*, int*, passivedouble*, passivedouble*, int*, - passivedouble*, int*, passivedouble*, passivedouble*, int*); -#endif - using namespace std; diff --git a/Common/include/linear_algebra/CMatrixVectorProduct.hpp b/Common/include/linear_algebra/CMatrixVectorProduct.hpp new file mode 100644 index 000000000000..a15c35f3926f --- /dev/null +++ b/Common/include/linear_algebra/CMatrixVectorProduct.hpp @@ -0,0 +1,167 @@ +/*! + * \file CMatrixVectorProduct.hpp + * \brief Headers for the classes related to sparse matrix-vector product wrappers. + * The actual operations are currently implemented mostly by CSysMatrix. + * \author F. Palacios, J. Hicken, T. Economon + * \version 6.2.0 "Falcon" + * + * The current SU2 release has been coordinated by the + * SU2 International Developers Society + * with selected contributions from the open-source community. + * + * The main research teams contributing to the current release are: + * - Prof. Juan J. Alonso's group at Stanford University. + * - Prof. Piero Colonna's group at Delft University of Technology. + * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. + * - Prof. Alberto Guardone's group at Polytechnic University of Milan. + * - Prof. Rafael Palacios' group at Imperial College London. + * - Prof. Vincent Terrapon's group at the University of Liege. + * - Prof. Edwin van der Weide's group at the University of Twente. + * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. + * + * Copyright 2012-2019, Francisco D. Palacios, Thomas D. Economon, + * Tim Albring, and the SU2 contributors. + * + * 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 . + */ + +#pragma once + +#include "../config_structure.hpp" +#include "../vector_structure.hpp" +#include "../matrix_structure.hpp" +#include "../geometry_structure.hpp" + +/*! + * \class CMatrixVectorProduct + * \brief Abstract base class for defining matrix-vector products + * \author J. Hicken. + * + * The Krylov-subspace solvers require only matrix-vector products and + * not the actual matrix/Jacobian. We need some way to indicate which + * function will perform the product. However, sometimes the + * functions that define the product will require different numbers + * and types of inputs. For example, the forward-difference + * approximation to a Jacobian-vector product requires the vector that + * defines the Jacobian and a perturbation parameter. The + * CMatrixVectorProduct class is used to derive child classes that can + * handle the different types of matrix-vector products and still be + * passed to a single implementation of the Krylov solvers. + * This abstraction may also be used to define matrix-free products. + */ +template +class CMatrixVectorProduct { +public: + virtual ~CMatrixVectorProduct() = 0; ///< class destructor + virtual void operator()(const CSysVector & u, CSysVector & v) + const = 0; ///< matrix-vector product operation +}; +template +CMatrixVectorProduct::~CMatrixVectorProduct() {} + + +/*! + * \class CSysMatrixVectorProduct + * \brief Specialization of matrix-vector product that uses CSysMatrix class + */ +template +class CSysMatrixVectorProduct : public CMatrixVectorProduct { +private: + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the product. */ + CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ + CConfig* config; /*!< \brief pointer to matrix that defines the config. */ + + /*! + * \brief Default constructor of the class + * \note This class cannot be default constructed as that would leave us with invalid pointers. + */ + CSysMatrixVectorProduct(); + +public: + + /*! + * \brief constructor of the class + * \param[in] matrix_ref - matrix reference that will be used to define the products + * \param[in] geometry_ref - geometry associated with the problem + * \param[in] config_ref - config of the problem + */ + inline CSysMatrixVectorProduct(CSysMatrix & matrix_ref, + CGeometry *geometry_ref, CConfig *config_ref) { + sparse_matrix = &matrix_ref; + geometry = geometry_ref; + config = config_ref; + } + + /*! + * \brief destructor of the class + */ + ~CSysMatrixVectorProduct() {} + + /*! + * \brief operator that defines the CSysMatrix-CSysVector product + * \param[in] u - CSysVector that is being multiplied by the sparse matrix + * \param[out] v - CSysVector that is the result of the product + */ + inline void operator()(const CSysVector & u, CSysVector & v) const { + sparse_matrix->MatrixVectorProduct(u, v, geometry, config); + } +}; + + +/*! + * \class CSysMatrixVectorProductTransposed + * \brief Specialization of matrix-vector product that uses CSysMatrix class for transposed products + */ +template +class CSysMatrixVectorProductTransposed : public CMatrixVectorProduct { +private: + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the product. */ + CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ + CConfig* config; /*!< \brief pointer to matrix that defines the config. */ + + /*! + * \brief Default constructor of the class + * \note This class cannot be default constructed as that would leave us with invalid pointers. + */ + CSysMatrixVectorProductTransposed(); + +public: + + /*! + * \brief constructor of the class + * \param[in] matrix_ref - matrix reference that will be used to define the products + * \param[in] geometry_ref - geometry associated with the problem + * \param[in] config_ref - config of the problem + */ + inline CSysMatrixVectorProductTransposed(CSysMatrix & matrix_ref, + CGeometry *geometry_ref, CConfig *config_ref) { + sparse_matrix = &matrix_ref; + geometry = geometry_ref; + config = config_ref; + } + + /*! + * \brief destructor of the class + */ + ~CSysMatrixVectorProductTransposed() {} + + /*! + * \brief operator that defines the CSysMatrix-CSysVector product + * \param[in] u - CSysVector that is being multiplied by the sparse matrix + * \param[out] v - CSysVector that is the result of the product + */ + inline void operator()(const CSysVector & u, CSysVector & v) const { + sparse_matrix->MatrixVectorProductTransposed(u, v, geometry, config); + } +}; diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp new file mode 100644 index 000000000000..308bbd8ca965 --- /dev/null +++ b/Common/include/linear_algebra/CPreconditioner.hpp @@ -0,0 +1,254 @@ +/*! + * \file CPreconditioner.hpp + * \brief Headers for the classes related to linear preconditioner wrappers. + * The actual operations are currently implemented mostly by CSysMatrix. + * \author F. Palacios, J. Hicken, T. Economon + * \version 6.2.0 "Falcon" + * + * The current SU2 release has been coordinated by the + * SU2 International Developers Society + * with selected contributions from the open-source community. + * + * The main research teams contributing to the current release are: + * - Prof. Juan J. Alonso's group at Stanford University. + * - Prof. Piero Colonna's group at Delft University of Technology. + * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. + * - Prof. Alberto Guardone's group at Polytechnic University of Milan. + * - Prof. Rafael Palacios' group at Imperial College London. + * - Prof. Vincent Terrapon's group at the University of Liege. + * - Prof. Edwin van der Weide's group at the University of Twente. + * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. + * + * Copyright 2012-2019, Francisco D. Palacios, Thomas D. Economon, + * Tim Albring, and the SU2 contributors. + * + * 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 . + */ + +#pragma once + +#include "../config_structure.hpp" +#include "../vector_structure.hpp" +#include "../matrix_structure.hpp" +#include "../geometry_structure.hpp" + +/*! + * \class CPreconditioner + * \brief abstract base class for defining preconditioning operation + * \author J. Hicken. + * + * See the remarks regarding the CMatrixVectorProduct class. The same + * idea applies here to the preconditioning operation. + */ +template +class CPreconditioner { +public: + virtual ~CPreconditioner() = 0; ///< class destructor + virtual void operator()(const CSysVector & u, CSysVector & v) + const = 0; ///< preconditioning operation +}; +template +CPreconditioner::~CPreconditioner() {} + + +/*! + * \class CJacobiPreconditioner + * \brief specialization of preconditioner that uses CSysMatrix class + */ +template +class CJacobiPreconditioner : public CPreconditioner { +private: + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ + CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ + CConfig* config; /*!< \brief pointer to matrix that defines the config. */ + + /*! + * \brief Default constructor of the class + * \note This class cannot be default constructed as that would leave us with invalid pointers. + */ + CJacobiPreconditioner(); + +public: + + /*! + * \brief constructor of the class + * \param[in] matrix_ref - matrix reference that will be used to define the preconditioner + * \param[in] geometry_ref - geometry associated with the problem + * \param[in] config_ref - config of the problem + */ + inline CJacobiPreconditioner(CSysMatrix & matrix_ref, + CGeometry *geometry_ref, CConfig *config_ref) { + sparse_matrix = &matrix_ref; + geometry = geometry_ref; + config = config_ref; + } + + /*! + * \brief destructor of the class + */ + ~CJacobiPreconditioner() {} + + /*! + * \brief operator that defines the preconditioner operation + * \param[in] u - CSysVector that is being preconditioned + * \param[out] v - CSysVector that is the result of the preconditioning + */ + inline void operator()(const CSysVector & u, CSysVector & v) const { + sparse_matrix->ComputeJacobiPreconditioner(u, v, geometry, config); + } +}; + + +/*! + * \class CILUPreconditioner + * \brief specialization of preconditioner that uses CSysMatrix class + */ +template +class CILUPreconditioner : public CPreconditioner { +private: + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ + CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ + CConfig* config; /*!< \brief pointer to matrix that defines the config. */ + + /*! + * \brief Default constructor of the class + * \note This class cannot be default constructed as that would leave us with invalid pointers. + */ + CILUPreconditioner(); + +public: + + /*! + * \brief constructor of the class + * \param[in] matrix_ref - matrix reference that will be used to define the preconditioner + * \param[in] geometry_ref - geometry associated with the problem + * \param[in] config_ref - config of the problem + */ + inline CILUPreconditioner(CSysMatrix & matrix_ref, + CGeometry *geometry_ref, CConfig *config_ref) { + sparse_matrix = &matrix_ref; + geometry = geometry_ref; + config = config_ref; + } + + /*! + * \brief destructor of the class + */ + ~CILUPreconditioner() {} + + /*! + * \brief operator that defines the preconditioner operation + * \param[in] u - CSysVector that is being preconditioned + * \param[out] v - CSysVector that is the result of the preconditioning + */ + inline void operator()(const CSysVector & u, CSysVector & v) const { + sparse_matrix->ComputeILUPreconditioner(u, v, geometry, config); + } +}; + + +/*! + * \class CLU_SGSPreconditioner + * \brief specialization of preconditioner that uses CSysMatrix class + */ +template +class CLU_SGSPreconditioner : public CPreconditioner { +private: + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ + CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ + CConfig* config; /*!< \brief pointer to matrix that defines the config. */ + + /*! + * \brief Default constructor of the class + * \note This class cannot be default constructed as that would leave us with invalid pointers. + */ + CLU_SGSPreconditioner(); + +public: + + /*! + * \brief constructor of the class + * \param[in] matrix_ref - matrix reference that will be used to define the preconditioner + * \param[in] geometry_ref - geometry associated with the problem + * \param[in] config_ref - config of the problem + */ + inline CLU_SGSPreconditioner(CSysMatrix & matrix_ref, + CGeometry *geometry_ref, CConfig *config_ref) { + sparse_matrix = &matrix_ref; + geometry = geometry_ref; + config = config_ref; + } + + /*! + * \brief destructor of the class + */ + ~CLU_SGSPreconditioner() {} + + /*! + * \brief operator that defines the preconditioner operation + * \param[in] u - CSysVector that is being preconditioned + * \param[out] v - CSysVector that is the result of the preconditioning + */ + inline void operator()(const CSysVector & u, CSysVector & v) const { + sparse_matrix->ComputeLU_SGSPreconditioner(u, v, geometry, config); + } +}; + + +/*! + * \class CLineletPreconditioner + * \brief specialization of preconditioner that uses CSysMatrix class + */ +template +class CLineletPreconditioner : public CPreconditioner { +private: + CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ + CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ + CConfig* config; /*!< \brief pointer to matrix that defines the config. */ + + /*! + * \brief Default constructor of the class + * \note This class cannot be default constructed as that would leave us with invalid pointers. + */ + CLineletPreconditioner(); + +public: + + /*! + * \brief constructor of the class + * \param[in] matrix_ref - matrix reference that will be used to define the preconditioner + * \param[in] geometry_ref - geometry associated with the problem + * \param[in] config_ref - config of the problem + */ + inline CLineletPreconditioner(CSysMatrix & matrix_ref, + CGeometry *geometry_ref, CConfig *config_ref) { + sparse_matrix = &matrix_ref; + geometry = geometry_ref; + config = config_ref; + } + + /*! + * \brief destructor of the class + */ + ~CLineletPreconditioner() {} + + /*! + * \brief operator that defines the preconditioner operation + * \param[in] u - CSysVector that is being preconditioned + * \param[out] v - CSysVector that is the result of the preconditioning + */ + inline void operator()(const CSysVector & u, CSysVector & v) const { + sparse_matrix->ComputeLineletPreconditioner(u, v, geometry, config); + } +}; diff --git a/Common/include/linear_solvers_structure.hpp b/Common/include/linear_solvers_structure.hpp index 2133a6bf26c0..663785db2ccc 100644 --- a/Common/include/linear_solvers_structure.hpp +++ b/Common/include/linear_solvers_structure.hpp @@ -54,6 +54,8 @@ #include "matrix_structure.hpp" #include "config_structure.hpp" #include "geometry_structure.hpp" +#include "linear_algebra/CMatrixVectorProduct.hpp" +#include "linear_algebra/CPreconditioner.hpp" using namespace std; @@ -85,6 +87,7 @@ class CSysSolve { bool cg_ready; /*!< \brief Indicate if memory used by CG is allocated. */ bool bcg_ready; /*!< \brief Indicate if memory used by BCGSTAB is allocated. */ bool gmres_ready; /*!< \brief Indicate if memory used by FGMRES is allocated. */ + bool smooth_ready; /*!< \brief Indicate if memory used by SMOOTHER is allocated. */ VectorType r; /*!< \brief Residual in CG and BCGSTAB. */ VectorType A_x; /*!< \brief Result of matrix-vector product in CG and BCGSTAB. */ @@ -111,7 +114,10 @@ class CSysSolve { * so, feel free to delete this and replace it as needed with the * appropriate global function */ - ScalarType Sign(const ScalarType & x, const ScalarType & y) const; + inline ScalarType Sign(const ScalarType & x, const ScalarType & y) const { + if (y == 0.0) return 0.0; + return fabs(x) * (y < 0.0 ? -1.0 : 1.0); + } /*! * \brief applys a Givens rotation to a 2-vector @@ -211,7 +217,7 @@ class CSysSolve { /*! \brief Conjugate Gradient method * \param[in] b - the right hand size vector - * \param[in, out] x - on entry the intial guess, on exit the solution + * \param[in,out] x - on entry the intial guess, on exit the solution * \param[in] mat_vec - object that defines matrix-vector product * \param[in] precond - object that defines preconditioner * \param[in] tol - tolerance with which to solve the system @@ -226,12 +232,12 @@ class CSysSolve { /*! * \brief Flexible Generalized Minimal Residual method * \param[in] b - the right hand size vector - * \param[in, out] x - on entry the intial guess, on exit the solution + * \param[in,out] x - on entry the intial guess, on exit the solution * \param[in] mat_vec - object that defines matrix-vector product * \param[in] precond - object that defines preconditioner * \param[in] tol - tolerance with which to solve the system * \param[in] m - maximum size of the search subspace - * \param[in] residual + * \param[in] residual - norm of final residual * \param[in] monitoring - turn on priting residuals from solver to screen. * \param[in] config - Definition of the particular problem. */ @@ -242,12 +248,12 @@ class CSysSolve { /*! * \brief Biconjugate Gradient Stabilized Method (BCGSTAB) * \param[in] b - the right hand size vector - * \param[in, out] x - on entry the intial guess, on exit the solution + * \param[in,out] x - on entry the intial guess, on exit the solution * \param[in] mat_vec - object that defines matrix-vector product * \param[in] precond - object that defines preconditioner * \param[in] tol - tolerance with which to solve the system * \param[in] m - maximum size of the search subspace - * \param[in] residual + * \param[in] residual - norm of final residual * \param[in] monitoring - turn on priting residuals from solver to screen. * \param[in] config - Definition of the particular problem. */ @@ -255,6 +261,22 @@ class CSysSolve { PrecondType & precond, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CConfig *config); + /*! + * \brief Generic smoother (modified Richardson iteration with preconditioner) + * \param[in] b - the right hand size vector + * \param[in,out] x - on entry the intial guess, on exit the solution + * \param[in] mat_vec - object that defines matrix-vector product + * \param[in] precond - object that defines preconditioner + * \param[in] tol - tolerance with which to solve the system + * \param[in] m - maximum number of iterations + * \param[in] residual - norm of final residual + * \param[in] monitoring - turn on priting residuals from solver to screen. + * \param[in] config - Definition of the particular problem. + */ + unsigned long Smoother_LinSolver(const VectorType & b, VectorType & x, ProductType & mat_vec, + PrecondType & precond, ScalarType tol, unsigned long m, + ScalarType *residual, bool monitoring, CConfig *config); + /*! * \brief Solve the linear system using a Krylov subspace method * \param[in] Jacobian - Jacobian Matrix for the linear system @@ -281,8 +303,6 @@ class CSysSolve { * \brief Get the final residual. * \return The residual at the end of Solve */ - ScalarType GetResidual(void) const; + inline ScalarType GetResidual(void) const { return Residual; } }; - -#include "linear_solvers_structure.inl" diff --git a/Common/include/linear_solvers_structure.inl b/Common/include/linear_solvers_structure.inl deleted file mode 100644 index b229828027df..000000000000 --- a/Common/include/linear_solvers_structure.inl +++ /dev/null @@ -1,58 +0,0 @@ -/*! - * \file linear_solvers_structure.inl - * \brief inline subroutines of the linear_solvers_structure.hpp file. - * \author J. Hicken, F. Palacios, T. Economon - * \version 6.2.0 "Falcon" - * - * The current SU2 release has been coordinated by the - * SU2 International Developers Society - * with selected contributions from the open-source community. - * - * The main research teams contributing to the current release are: - * - Prof. Juan J. Alonso's group at Stanford University. - * - Prof. Piero Colonna's group at Delft University of Technology. - * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. - * - Prof. Alberto Guardone's group at Polytechnic University of Milan. - * - Prof. Rafael Palacios' group at Imperial College London. - * - Prof. Vincent Terrapon's group at the University of Liege. - * - Prof. Edwin van der Weide's group at the University of Twente. - * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. - * - * Copyright 2012-2019, Francisco D. Palacios, Thomas D. Economon, - * Tim Albring, and the SU2 contributors. - * - * 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 . - */ - -#pragma once - -template -inline ScalarType CSysSolve::Sign(const ScalarType & x, const ScalarType & y) const { - if (y == 0.0) - return 0.0; - else { -// return (y < 0 ? -fabs(x) : fabs(x)); - if (y < 0) return -fabs(x); - else return fabs(x); - } -} - -template -inline ScalarType CSysSolve::GetResidual(void) const { return Residual; } - -template -void CSysSolve::HandleTemporariesIn(CSysVector & LinSysRes, CSysVector & LinSysSol) {} - -template -void CSysSolve::HandleTemporariesOut(CSysVector & LinSysSol) {} diff --git a/Common/include/matrix_structure.hpp b/Common/include/matrix_structure.hpp index c6db8c5be5a7..72f15266ca49 100644 --- a/Common/include/matrix_structure.hpp +++ b/Common/include/matrix_structure.hpp @@ -48,8 +48,16 @@ #include "geometry_structure.hpp" #include "vector_structure.hpp" -#ifdef HAVE_MKL +#if defined(HAVE_MKL) && !defined(CODI_FORWARD_TYPE) #include "mkl.h" +#define USE_MKL +/*--- + Lapack direct calls only seem to be created for Intel compilers, and it is not worthwhile + making "getrf" and "getrs" compatible with AD since they are not used as often as "gemm". +---*/ +#if defined(__INTEL_COMPILER) && defined(MKL_DIRECT_CALL_SEQ) && !defined(CODI_REVERSE_TYPE) + #define USE_MKL_LAPACK +#endif #endif using namespace std; @@ -83,68 +91,55 @@ class CSysMatrix { unsigned short ilu_fill_in; /*!< \brief Fill in level for the ILU preconditioner. */ ScalarType *block; /*!< \brief Internal array to store a subblock of the matrix. */ - ScalarType *block_inverse; /*!< \brief Internal array to store a subblock of the matrix. */ - ScalarType *block_weight; /*!< \brief Internal array to store a subblock of the matrix. */ - ScalarType *prod_block_vector; /*!< \brief Internal array to store the product of a subblock with a vector. */ + ScalarType *block_inverse; /*!< \brief Internal array to store a subblock of the matrix. */ + ScalarType *block_weight; /*!< \brief Internal array to store a subblock of the matrix. */ ScalarType *prod_row_vector; /*!< \brief Internal array to store the product of a matrix-by-blocks "row" with a vector. */ - ScalarType *aux_vector; /*!< \brief Auxiliary array to store intermediate results. */ - ScalarType *sum_vector; /*!< \brief Auxiliary array to store intermediate results. */ - ScalarType *invM; /*!< \brief Inverse of (Jacobi) preconditioner. */ - - bool *LineletBool; /*!< \brief Identify if a point belong to a linelet. */ - vector *LineletPoint; /*!< \brief Linelet structure. */ - unsigned long nLinelet; /*!< \brief Number of Linelets in the system. */ - ScalarType **UBlock, **invUBlock, **LBlock, - **yVector, **zVector, **rVector, *LFBlock, - *LyVector, *FzVector; /*!< \brief Arrays of the Linelet preconditioner methodology. */ - unsigned long max_nElem; + ScalarType *aux_vector; /*!< \brief Auxiliary array to store intermediate results. */ + ScalarType *sum_vector; /*!< \brief Auxiliary array to store intermediate results. */ + ScalarType *invM; /*!< \brief Inverse of (Jacobi) preconditioner, or diagonal of ILU. */ + + unsigned long nLinelet; /*!< \brief Number of Linelets in the system. */ + vector LineletBool; /*!< \brief Identify if a point belong to a Linelet. */ + vector > LineletPoint; /*!< \brief Linelet structure. */ + vector LineletUpper; /*!< \brief Pointers to the upper blocks of the tri-diag system. */ + vector LineletInvDiag; /*!< \brief Inverse of the diagonal blocks of the tri-diag system. */ + vector LineletVector; /*!< \brief Solution and RHS of the tri-diag system. */ -#if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) - void * MatrixMatrixProductJitter; /*!< \brief Jitter handle for MKL JIT based GEMM. */ - dgemm_jit_kernel_t MatrixMatrixProductKernel; /*!< \brief MKL JIT based GEMM kernel. */ - void * MatrixVectorProductJitterBetaZero; /*!< \brief Jitter handle for MKL JIT based GEMV. */ - dgemm_jit_kernel_t MatrixVectorProductKernelBetaZero; /*!< \brief MKL JIT based GEMV kernel. */ - void * MatrixVectorProductJitterBetaOne; /*!< \brief Jitter handle for MKL JIT based GEMV with BETA=1.0. */ - dgemm_jit_kernel_t MatrixVectorProductKernelBetaOne; /*!< \brief MKL JIT based GEMV kernel with BETA=1.0. */ - bool useMKL; +#ifdef USE_MKL + void * MatrixMatrixProductJitter; /*!< \brief Jitter handle for MKL JIT based GEMM. */ + dgemm_jit_kernel_t MatrixMatrixProductKernel; /*!< \brief MKL JIT based GEMM kernel. */ + void * MatrixVectorProductJitterBetaZero; /*!< \brief Jitter handle for MKL JIT based GEMV. */ + dgemm_jit_kernel_t MatrixVectorProductKernelBetaZero; /*!< \brief MKL JIT based GEMV kernel. */ + void * MatrixVectorProductJitterBetaOne; /*!< \brief Jitter handle for MKL JIT based GEMV with BETA=1.0. */ + dgemm_jit_kernel_t MatrixVectorProductKernelBetaOne; /*!< \brief MKL JIT based GEMV kernel with BETA=1.0. */ + void * MatrixVectorProductJitterAlphaMinusOne; /*!< \brief Jitter handle for MKL JIT based GEMV with ALPHA=-1.0 and BETA=1.0. */ + dgemm_jit_kernel_t MatrixVectorProductKernelAlphaMinusOne; /*!< \brief MKL JIT based GEMV kernel with ALPHA=-1.0 and BETA=1.0. */ + void * MatrixVectorProductTranspJitterBetaOne; /*!< \brief Jitter handle for MKL JIT based GEMV (transposed) with BETA=1.0. */ + dgemm_jit_kernel_t MatrixVectorProductTranspKernelBetaOne; /*!< \brief MKL JIT based GEMV (transposed) kernel with BETA=1.0. */ + lapack_int * mkl_ipiv; #endif /*! * \brief Handle type conversion for when we Set, Add, etc. blocks, preserving derivative information (if supported by types). + * \note See specializations for discrete adjoint right outside this class's declaration. */ template - DstType ActiveAssign(const SrcType & val) const; + inline DstType ActiveAssign(const SrcType & val) const { return val; } /*! * \brief Handle type conversion for when we Set, Add, etc. blocks, discarding derivative information. */ template - DstType PassiveAssign(const SrcType & val) const; - -public: - - /*! - * \brief Constructor of the class. - */ - CSysMatrix(void); - - /*! - * \brief Destructor of the class. - */ - ~CSysMatrix(void); - - /*! - * \brief Initializes space matrix system. - * \param[in] nVar - Number of variables. - * \param[in] nEqn - Number of equations. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void Initialize(unsigned long nPoint, unsigned long nPointDomain, unsigned short nVar, unsigned short nEqn, - bool EdgeConnect, CGeometry *geometry, CConfig *config); + inline DstType PassiveAssign(const SrcType & val) const { +#if defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE) + return SU2_TYPE::GetValue(val); +#else + return val; +#endif + } /*! - * \brief Assigns values to the sparse-matrix structure. + * \brief Assigns values to the sparse-matrix structure (used in Initialize). * \param[in] val_nPoint - Number of points in the nPoint x nPoint block structure * \param[in] val_nVar - Number of nVar x nVar variables in each subblock of the matrix-by-block structure. * \param[in] val_nEq - Number of nEqn x nVar variables in each subblock of the matrix-by-block structure. @@ -156,7 +151,7 @@ class CSysMatrix { void SetIndexes(unsigned long val_nPoint, unsigned long val_nPointDomain, unsigned short val_nVar, unsigned short val_nEq, unsigned long* val_row_ptr, unsigned long* val_col_ind, unsigned long val_nnz, CConfig *config); /*! - * \brief Assigns values to the sparse-matrix structure. + * \brief Assigns values to the sparse-matrix structure (used in Initialize). * \param[in] geometry - Geometrical definition of the problem. * \param[in] iPoint - Base point to compute neighbours. * \param[in] deep_level - Deep level for the recursive algorithm. @@ -167,92 +162,104 @@ class CSysMatrix { void SetNeighbours(CGeometry *geometry, unsigned long iPoint, unsigned short deep_level, unsigned short fill_level, bool EdgeConnect, vector & vneighs); /*! - * \brief Sets to zero all the entries of the sparse matrix. + * \brief Calculates the matrix-vector product: product = matrix*vector + * \param[in] matrix + * \param[in] vector + * \param[out] product */ - void SetValZero(void); + inline void MatrixVectorProduct(const ScalarType *matrix, const ScalarType *vector, ScalarType *product); /*! - * \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. + * \brief Calculates the matrix-vector product: product += matrix*vector + * \param[in] matrix + * \param[in] vector + * \param[in,out] product */ - template - void InitiateComms(CSysVector & x, - CGeometry *geometry, - CConfig *config, - unsigned short commType); + inline void MatrixVectorProductAdd(const ScalarType *matrix, const ScalarType *vector, ScalarType *product); /*! - * \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. + * \brief Calculates the matrix-vector product: product -= matrix*vector + * \param[in] matrix + * \param[in] vector + * \param[in,out] product */ - template - void CompleteComms(CSysVector & x, - CGeometry *geometry, - CConfig *config, - unsigned short commType); + inline void MatrixVectorProductSub(const ScalarType *matrix, const ScalarType *vector, ScalarType *product); /*! - * \brief Copies the block (i, j) of the matrix-by-blocks structure in the internal variable *block. - * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. - * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. + * \brief Calculates the matrix-vector product: product += matrix^T * vector + * \param[in] matrix + * \param[in] vector + * \param[in,out] product */ - ScalarType *GetBlock(unsigned long block_i, unsigned long block_j); + inline void MatrixVectorProductTransp(const ScalarType *matrix, const ScalarType *vector, ScalarType *product); /*! - * \brief Copies the block (i, j) of the matrix-by-blocks structure in the internal variable *block. - * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. - * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. + * \brief Calculates the matrix-matrix product + * \param[in] matrix_a + * \param[in] matrix_b + * \param[out] product */ - ScalarType GetBlock(unsigned long block_i, unsigned long block_j, unsigned short iVar, unsigned short jVar); + inline void MatrixMatrixProduct(const ScalarType *matrix_a, const ScalarType *matrix_b, ScalarType *product); /*! - * \brief Set the value of a block in the sparse matrix. - * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. - * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. - * \param[in] **val_block - Block to set to A(i, j). + * \brief Subtract b from a and store the result in c. */ - template - void SetBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block); + inline void VectorSubtraction(const ScalarType *a, const ScalarType *b, ScalarType *c) { + for(unsigned long iVar = 0; iVar < nVar; iVar++) + c[iVar] = a[iVar] - b[iVar]; + } /*! - * \brief Set the value of a block in the sparse matrix. - * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. - * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. - * \param[in] **val_block - Block to set to A(i, j). + * \brief Subtract b from a and store the result in c. */ - template - void SetBlock(unsigned long block_i, unsigned long block_j, OtherType *val_block); + inline void MatrixSubtraction(const ScalarType *a, const ScalarType *b, ScalarType *c) { + for(unsigned long iVar = 0; iVar < nVar*nEqn; iVar++) + c[iVar] = a[iVar] - b[iVar]; + } /*! - * \brief Adds the specified block to the sparse matrix. + * \brief Solve a small (nVar x nVar) linear system using Gaussian elimination. + * \param[in,out] matrix - On entry the system matrix, on exit the factorized matrix. + * \param[in,out] vec - On entry the rhs, on exit the solution. + */ + inline void Gauss_Elimination(ScalarType* matrix, ScalarType* vec); + + /*! + * \brief Invert a small dense matrix. + * \param[in] matrix - the matrix. + * \param[out] inverse - the matrix inverse. + */ + inline void MatrixInverse(const ScalarType *matrix, ScalarType *inverse); + + /*! + * \brief Performs the Gauss Elimination algorithm to solve the linear subsystem of the (i, i) subblock and rhs. + * \param[in] block_i - Index of the (i, i) subblock in the matrix-by-blocks structure. + * \param[in] rhs - Right-hand-side of the linear system. + * \param[in] transposed - If true the transposed of the block is used (default = false). + * \return Solution of the linear system (overwritten on rhs). + */ + inline void Gauss_Elimination(unsigned long block_i, ScalarType* rhs, bool transposed = false); + + /*! + * \brief Inverse diagonal block. * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. - * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. - * \param[in] **val_block - Block to add to A(i, j). + * \param[out] invBlock - Inverse block. */ - template - void AddBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block); + inline void InverseDiagonalBlock(unsigned long block_i, ScalarType *invBlock, bool transpose = false); /*! - * \brief Subtracts the specified block to the sparse matrix. + * \brief Inverse diagonal block. * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. - * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. - * \param[in] **val_block - Block to subtract to A(i, j). + * \param[out] invBlock - Inverse block. */ - template - void SubtractBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block); + inline void InverseDiagonalBlock_ILUMatrix(unsigned long block_i, ScalarType *invBlock); /*! * \brief Copies the block (i, j) of the matrix-by-blocks structure in the internal variable *block. * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. */ - ScalarType *GetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j); + inline ScalarType *GetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j); /*! * \brief Set the value of a block in the sparse matrix. @@ -260,7 +267,7 @@ class CSysMatrix { * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. * \param[in] **val_block - Block to set to A(i, j). */ - void SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block); + inline void SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block); /*! @@ -269,7 +276,7 @@ class CSysMatrix { * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. * \param[in] **val_block - Block to set to A(i, j). */ - void SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block); + inline void SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block); /*! * \brief Subtracts the specified block to the sparse matrix. @@ -277,204 +284,291 @@ class CSysMatrix { * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. * \param[in] **val_block - Block to subtract to A(i, j). */ - void SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block); + inline void SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block); /*! - * \brief Adds the specified value to the diagonal of the (i, i) subblock - * of the matrix-by-blocks structure. - * \param[in] block_i - Index of the block in the matrix-by-blocks structure. - * \param[in] val_matrix - Value to add to the diagonal elements of A(i, i). - */ - template - void AddVal2Diag(unsigned long block_i, OtherType val_matrix); - - /*! - * \brief Sets the specified value to the diagonal of the (i, i) subblock - * of the matrix-by-blocks structure. - * \param[in] block_i - Index of the block in the matrix-by-blocks structure. - * \param[in] val_matrix - Value to add to the diagonal elements of A(i, i). + * \brief Performs the product of i-th row of the upper part of a sparse matrix by a vector. + * \param[in] vec - Vector to be multiplied by the upper part of the sparse matrix A. + * \param[in] row_i - Row of the matrix to be multiplied by vector vec. + * \return prod Result of the product U(A)*vec (stored at *prod_row_vector). */ - template - void SetVal2Diag(unsigned long block_i, OtherType val_matrix); + void UpperProduct(const CSysVector & vec, unsigned long row_i); /*! - * \brief Calculates the matrix-vector product - * \param[in] matrix - * \param[in] vector - * \param[out] product + * \brief Performs the product of i-th row of the lower part of a sparse matrix by a vector. + * \param[in] vec - Vector to be multiplied by the lower part of the sparse matrix A. + * \param[in] row_i - Row of the matrix to be multiplied by vector vec. + * \return prod Result of the product L(A)*vec (stored at *prod_row_vector). */ - void MatrixVectorProduct(ScalarType *matrix, ScalarType *vector, ScalarType *product); + void LowerProduct(const CSysVector & vec, unsigned long row_i); /*! - * \brief Calculates the matrix-matrix product - * \param[in] matrix_a - * \param[in] matrix_b - * \param[out] product + * \brief Performs the product of i-th row of the diagonal part of a sparse matrix by a vector. + * \param[in] vec - Vector to be multiplied by the diagonal part of the sparse matrix A. + * \param[in] row_i - Row of the matrix to be multiplied by vector vec. + * \return prod Result of the product D(A)*vec (stored at *prod_row_vector). */ - void MatrixMatrixProduct(ScalarType *matrix_a, ScalarType *matrix_b, ScalarType *product); + void DiagonalProduct(const CSysVector & vec, unsigned long row_i); /*! - * \brief Deletes the values of the row i of the sparse matrix. - * \param[in] i - Index of the row. + * \brief Performs the product of i-th row of a sparse matrix by a vector. + * \param[in] vec - Vector to be multiplied by the row of the sparse matrix A. + * \param[in] row_i - Row of the matrix to be multiplied by vector vec. + * \return Result of the product (stored at *prod_row_vector). */ - void DeleteValsRowi(unsigned long i); + void RowProduct(const CSysVector & vec, unsigned long row_i); - /*! - * \brief Recursive definition of determinate using expansion by minors. Written by Paul Bourke - * \param[in] a - Matrix to compute the determinant. - * \param[in] n - Size of the quare matrix. - * \return Value of the determinant. - */ - ScalarType MatrixDeterminant(ScalarType **a, unsigned long n); +public: /*! - * \brief Find the cofactor matrix of a square matrix. Written by Paul Bourke - * \param[in] a - Matrix to compute the determinant. - * \param[in] n - Size of the quare matrix. - * \param[out] b - cofactor matrix + * \brief Constructor of the class. */ - void MatrixCoFactor(ScalarType **a, unsigned long n, ScalarType **b) ; + CSysMatrix(void); /*! - * \brief Transpose of a square matrix, do it in place. Written by Paul Bourke - * \param[in] a - Matrix to compute the determinant. - * \param[in] n - Size of the quare matrix. + * \brief Destructor of the class. */ - void MatrixTranspose(ScalarType **a, unsigned long n) ; + ~CSysMatrix(void); /*! - * \brief Performs the Gauss Elimination algorithm to solve the linear subsystem of the (i, i) subblock and rhs. - * \param[in] block_i - Index of the (i, i) subblock in the matrix-by-blocks structure. - * \param[in] rhs - Right-hand-side of the linear system. - * \param[in] transposed - If true the transposed of the block is used (default = false). - * \return Solution of the linear system (overwritten on rhs). + * \brief Initializes sparse matrix system. + * \param[in] nVar - Number of variables. + * \param[in] nEqn - Number of equations. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. */ - void Gauss_Elimination(unsigned long block_i, ScalarType* rhs, bool transposed = false); + void Initialize(unsigned long nPoint, unsigned long nPointDomain, unsigned short nVar, unsigned short nEqn, + bool EdgeConnect, CGeometry *geometry, CConfig *config); /*! - * \brief Performs the Gauss Elimination algorithm to solve the linear subsystem of the (i, i) subblock and rhs. - * \param[in] Block - matrix-by-blocks structure. - * \param[in] rhs - Right-hand-side of the linear system. - * \return Solution of the linear system (overwritten on rhs). + * \brief Sets to zero all the entries of the sparse matrix. */ - void Gauss_Elimination(ScalarType* Block, ScalarType* rhs); + inline void SetValZero(void) { + if(matrix != NULL) + for (unsigned long index = 0; index < nnz*nVar*nEqn; index++) + matrix[index] = 0.0; + } /*! - * \brief Performs the Gauss Elimination algorithm to solve the linear subsystem of the (i, i) subblock and rhs. - * \param[in] block_i - Index of the (i, i) subblock in the matrix-by-blocks structure. - * \param[in] rhs - Right-hand-side of the linear system. - * \return Solution of the linear system (overwritten on rhs). + * \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. */ - void Gauss_Elimination_ILUMatrix(unsigned long block_i, ScalarType* rhs); + template + void InitiateComms(CSysVector & x, + CGeometry *geometry, + CConfig *config, + unsigned short commType); /*! - * \fn void CSysMatrix::ProdBlockVector(unsigned long block_i, unsigned long block_j, su2double* vec); - * \brief Performs the product of the block (i, j) by vector vec. - * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. - * \param[in] block_j - Indexes of the block in the matrix-by-blocks structure. - * \param[in] vec - Vector to be multiplied by the block (i, j) of the sparse matrix A. - * \return Product of A(i, j) by vector *vec (stored at *prod_block_vector). + * \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. */ - void ProdBlockVector(unsigned long block_i, unsigned long block_j, const CSysVector & vec); + template + void CompleteComms(CSysVector & x, + CGeometry *geometry, + CConfig *config, + unsigned short commType); /*! - * \brief Performs the product of i-th row of the upper part of a sparse matrix by a vector. - * \param[in] vec - Vector to be multiplied by the upper part of the sparse matrix A. - * \param[in] row_i - Row of the matrix to be multiplied by vector vec. - * \return prod Result of the product U(A)*vec (stored at *prod_row_vector). + * \brief Get a pointer to the start of block "ij" + * \param[in] block_i - Row index. + * \param[in] block_j - Column index. + * \return Pointer to location in memory where the block starts. */ - void UpperProduct(CSysVector & vec, unsigned long row_i); + inline ScalarType *GetBlock(unsigned long block_i, unsigned long block_j) { + + for (unsigned long index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) + if (col_ind[index] == block_j) + return &(matrix[index*nVar*nEqn]); + + return NULL; + } /*! - * \brief Performs the product of i-th row of the lower part of a sparse matrix by a vector. - * \param[in] vec - Vector to be multiplied by the lower part of the sparse matrix A. - * \param[in] row_i - Row of the matrix to be multiplied by vector vec. - * \return prod Result of the product L(A)*vec (stored at *prod_row_vector). + * \brief Gets the value of a particular entry in block "ij". + * \param[in] block_i - Row index. + * \param[in] block_j - Column index. + * \param[in] iVar - Row of the block. + * \param[in] jVar - Column of the block. + * \return Value of the block entry. */ - void LowerProduct(CSysVector & vec, unsigned long row_i); + inline ScalarType GetBlock(unsigned long block_i, unsigned long block_j, + unsigned short iVar, unsigned short jVar) { - /*! - * \brief Performs the product of i-th row of the diagonal part of a sparse matrix by a vector. - * \param[in] vec - Vector to be multiplied by the diagonal part of the sparse matrix A. - * \param[in] row_i - Row of the matrix to be multiplied by vector vec. - * \return prod Result of the product D(A)*vec (stored at *prod_row_vector). - */ - void DiagonalProduct(CSysVector & vec, unsigned long row_i); + for (unsigned long index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) + if (col_ind[index] == block_j) + return matrix[index*nVar*nEqn+iVar*nEqn+jVar]; + + return 0.0; + } /*! - * \brief Performs the product of i-th row of a sparse matrix by a vector. - * \param[in] vec - Vector to be multiplied by the row of the sparse matrix A. - * \param[in] row_i - Row of the matrix to be multiplied by vector vec. - * \return Result of the product (stored at *prod_row_vector). + * \brief Set the value of a block in the sparse matrix. + * \param[in] block_i - Row index. + * \param[in] block_j - Column index. + * \param[in] **val_block - Block to set to A(i, j). */ - void RowProduct(const CSysVector & vec, unsigned long row_i); + template + inline void SetBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block) { + + unsigned long iVar, jVar, index; - /*! - * \brief Performs the product of a sparse matrix by a vector. - * \param[in] vec - Vector to be multiplied by the sparse matrix A. - * \param[out] prod - Result of the product. - * \return Result of the product A*vec. - */ - void MatrixVectorProduct(const CSysVector & vec, CSysVector & prod); + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + if (col_ind[index] == block_j) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++) + matrix[index*nVar*nEqn+iVar*nEqn+jVar] = PassiveAssign(val_block[iVar][jVar]); + break; + } + } + } /*! - * \brief Performs the product of a sparse matrix by a CSysVector. - * \param[in] vec - CSysVector to be multiplied by the sparse matrix A. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[out] prod - Result of the product. + * \brief Set the value of a block in the sparse matrix. + * \param[in] block_i - Row index. + * \param[in] block_j - Column index. + * \param[in] *val_block - Block to set to A(i, j). */ - void MatrixVectorProduct(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); + template + inline void SetBlock(unsigned long block_i, unsigned long block_j, OtherType *val_block) { + + unsigned long iVar, jVar, index; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + if (col_ind[index] == block_j) { + for (iVar = 0; iVar < nVar*nEqn; iVar++) + matrix[index*nVar*nEqn+iVar] = PassiveAssign(val_block[iVar]); + break; + } + } + } /*! - * \brief Performs the product of a sparse matrix by a CSysVector. - * \param[in] vec - CSysVector to be multiplied by the sparse matrix A. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[out] prod - Result of the product. + * \brief Adds the specified block to the sparse matrix. + * \param[in] block_i - Row index. + * \param[in] block_j - Column index. + * \param[in] **val_block - Block to add to A(i, j). */ - void MatrixVectorProductTransposed(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); + template + inline void AddBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block) { + + unsigned long iVar, jVar, index; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + if (col_ind[index] == block_j) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++) + matrix[index*nVar*nEqn+iVar*nEqn+jVar] += PassiveAssign(val_block[iVar][jVar]); + break; + } + } + } /*! - * \brief Performs the product of two block matrices. + * \brief Subtracts the specified block to the sparse matrix. + * \param[in] block_i - Row index. + * \param[in] block_j - Column index. + * \param[in] **val_block - Block to subtract to A(i, j). */ - void GetMultBlockBlock(ScalarType *c, ScalarType *a, ScalarType *b); + template + inline void SubtractBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block) { + + unsigned long iVar, jVar, index; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + if (col_ind[index] == block_j) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++) + matrix[index*nVar*nEqn+iVar*nEqn+jVar] -= PassiveAssign(val_block[iVar][jVar]); + break; + } + } + } /*! - * \brief Performs the product of a block matrices by a vector. + * \brief Adds the specified value to the diagonal of the (i, i) subblock + * of the matrix-by-blocks structure. + * \param[in] block_i - Diagonal index. + * \param[in] val_matrix - Value to add to the diagonal elements of A(i, i). */ - void GetMultBlockVector(ScalarType *c, ScalarType *a, ScalarType *b); + template + inline void AddVal2Diag(unsigned long block_i, OtherType val_matrix) { + + unsigned long iVar, index; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + if (col_ind[index] == block_i) { // Only elements on the diagonal + for (iVar = 0; iVar < nVar; iVar++) + matrix[index*nVar*nVar+iVar*nVar+iVar] += PassiveAssign(val_matrix); + break; + } + } + } /*! - * \brief Performs the subtraction of two matrices. + * \brief Sets the specified value to the diagonal of the (i, i) subblock + * of the matrix-by-blocks structure. + * \param[in] block_i - Diagonal index. + * \param[in] val_matrix - Value to add to the diagonal elements of A(i, i). */ - void GetSubsBlock(ScalarType *c, ScalarType *a, ScalarType *b); + template + inline void SetVal2Diag(unsigned long block_i, OtherType val_matrix) { + + unsigned long iVar, jVar, index; + + for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { + if (col_ind[index] == block_i) { // Only elements on the diagonal + + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) + matrix[index*nVar*nVar+iVar*nVar+jVar] = 0.0; + + for (iVar = 0; iVar < nVar; iVar++) + matrix[index*nVar*nVar+iVar*nVar+iVar] = PassiveAssign(val_matrix); + + break; + } + } + } /*! - * \brief Performs the subtraction of two vectors. + * \brief Deletes the values of the row i of the sparse matrix. + * \param[in] i - Index of the row. */ - void GetSubsVector(ScalarType *c, ScalarType *a, ScalarType *b); + void DeleteValsRowi(unsigned long i); /*! - * \brief Inverse diagonal block. - * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. - * \param[out] invBlock - Inverse block. + * \brief Modifies this matrix (A) and a rhs vector (b) such that (A^-1 * b)_i = x_i. + * \param[in] node_i - Index of the node for which to enforce the solution of all DOF's. + * \param[in] x_i - Values to enforce (nVar sized). + * \param[in,out] b - The rhs vector (b := b - A_{*,i} * x_i; b_i = x_i). */ - void InverseDiagonalBlock(unsigned long block_i, ScalarType *invBlock, bool transpose = false); + template + void EnforceSolutionAtNode(const unsigned long node_i, const OtherType *x_i, CSysVector & b); - /*! - * \brief Inverse diagonal block. - * \param[in] block_i - Indexes of the block in the matrix-by-blocks structure. - * \param[out] invBlock - Inverse block. + /*! + * \brief Performs the product of a sparse matrix by a CSysVector. + * \param[in] vec - CSysVector to be multiplied by the sparse matrix A. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[out] prod - Result of the product. */ - void InverseDiagonalBlock_ILUMatrix(unsigned long block_i, ScalarType *invBlock); + void MatrixVectorProduct(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); /*! - * \brief Inverse a block. - * \param[in] Block - block matrix. - * \param[out] invBlock - Inverse block. + * \brief Performs the product of a sparse matrix by a CSysVector. + * \param[in] vec - CSysVector to be multiplied by the sparse matrix A. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[out] prod - Result of the product. */ - void InverseBlock(ScalarType *Block, ScalarType *invBlock); + void MatrixVectorProductTransposed(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); /*! * \brief Build the Jacobi preconditioner. @@ -490,21 +584,6 @@ class CSysMatrix { */ void ComputeJacobiPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); - /*! - * \brief Apply Jacobi as a classical iterative smoother - * \param[in] b - CSysVector containing the residual (b) - * \param[in] x - CSysVector containing the solution (x^k) - * \param[in] mat_vec - object that defines matrix-vector product - * \param[in] tol - tolerance with which to solve the system - * \param[in] m - maximum size of the search subspace - * \param[in] residual - * \param[in] monitoring - turn on priting residuals from solver to screen. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[out] x - CSysVector containing the result of the smoothing (x^k+1 = x^k + M^-1*(b - A*x^k). - */ - unsigned long Jacobi_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config); - /*! * \brief Build the ILU preconditioner. * \param[in] transposed - Flag to use the transposed matrix to construct the preconditioner. @@ -520,21 +599,6 @@ class CSysMatrix { */ void ComputeILUPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); - /*! - * \brief Apply ILU as a classical iterative smoother - * \param[in] b - CSysVector containing the residual (b) - * \param[in] x - CSysVector containing the solution (x^k) - * \param[in] mat_vec - object that defines matrix-vector product - * \param[in] tol - tolerance with which to solve the system - * \param[in] m - maximum size of the search subspace - * \param[in] residual - * \param[in] monitoring - turn on priting residuals from solver to screen. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[out] x - CSysVector containing the result of the smoothing (x^k+1 = x^k + M^-1*(b - A*x^k). - */ - unsigned long ILU_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config); - /*! * \brief Multiply CSysVector by the preconditioner * \param[in] vec - CSysVector to be multiplied by the preconditioner. @@ -542,21 +606,6 @@ class CSysMatrix { */ void ComputeLU_SGSPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config); - /*! - * \brief Apply LU_SGS as a classical iterative smoother - * \param[in] b - CSysVector containing the residual (b) - * \param[in] x - CSysVector containing the solution (x^k) - * \param[in] mat_vec - object that defines matrix-vector product - * \param[in] tol - tolerance with which to solve the system - * \param[in] m - maximum size of the search subspace - * \param[in] residual - * \param[in] monitoring - turn on priting residuals from solver to screen. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[out] x - CSysVector containing the result of the smoothing (x^k+1 = x^k + M^-1*(b - A*x^k). - */ - unsigned long LU_SGS_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config); - /*! * \brief Build the Linelet preconditioner. * \param[in] geometry - Geometrical definition of the problem. @@ -581,242 +630,10 @@ class CSysMatrix { }; -/*! - * \class CSysMatrixVectorProduct - * \brief specialization of matrix-vector product that uses CSysMatrix class - */ -template -class CSysMatrixVectorProduct : public CMatrixVectorProduct { -private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the product. */ - CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ - CConfig* config; /*!< \brief pointer to matrix that defines the config. */ - -public: - - /*! - * \brief constructor of the class - * \param[in] matrix_ref - matrix reference that will be used to define the products - * \param[in] geometry_ref - - * \param[in] config_ref - - */ - CSysMatrixVectorProduct(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); - - /*! - * \brief destructor of the class - */ - ~CSysMatrixVectorProduct() {} - - /*! - * \brief operator that defines the CSysMatrix-CSysVector product - * \param[in] u - CSysVector that is being multiplied by the sparse matrix - * \param[out] v - CSysVector that is the result of the product - */ - void operator()(const CSysVector & u, CSysVector & v) const; -}; +#ifdef CODI_REVERSE_TYPE +template<> template<> +inline passivedouble CSysMatrix::ActiveAssign(const su2double & val) const { return SU2_TYPE::GetValue(val); } -/*! - * \class CSysMatrixVectorProduct - * \brief specialization of matrix-vector product that uses CSysMatrix class - */ -template -class CSysMatrixVectorProductTransposed : public CMatrixVectorProduct { -private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the product. */ - CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ - CConfig* config; /*!< \brief pointer to matrix that defines the config. */ - -public: - - /*! - * \brief constructor of the class - * \param[in] matrix_ref - matrix reference that will be used to define the products - * \param[in] geometry_ref - - * \param[in] config_ref - - */ - CSysMatrixVectorProductTransposed(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); - - /*! - * \brief destructor of the class - */ - ~CSysMatrixVectorProductTransposed() {} - - /*! - * \brief operator that defines the CSysMatrix-CSysVector product - * \param[in] u - CSysVector that is being multiplied by the sparse matrix - * \param[out] v - CSysVector that is the result of the product - */ - void operator()(const CSysVector & u, CSysVector & v) const; -}; - -/*! - * \class CJacobiPreconditioner - * \brief specialization of preconditioner that uses CSysMatrix class - */ -template -class CJacobiPreconditioner : public CPreconditioner { -private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ - CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ - CConfig* config; /*!< \brief pointer to matrix that defines the config. */ - -public: - - /*! - * \brief constructor of the class - * \param[in] matrix_ref - matrix reference that will be used to define the preconditioner - * \param[in] geometry_ref - - * \param[in] config_ref - - */ - CJacobiPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); - - /*! - * \brief destructor of the class - */ - ~CJacobiPreconditioner() {} - - /*! - * \brief operator that defines the preconditioner operation - * \param[in] u - CSysVector that is being preconditioned - * \param[out] v - CSysVector that is the result of the preconditioning - */ - void operator()(const CSysVector & u, CSysVector & v) const; -}; - -/*! - * \class CJacobiTransposedPreconditioner - * \brief specialization of preconditioner that uses CSysMatrix class - */ -template -class CJacobiTransposedPreconditioner : public CPreconditioner { -private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ - CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ - CConfig* config; /*!< \brief pointer to matrix that defines the config. */ - -public: - - /*! - * \brief constructor of the class - * \param[in] matrix_ref - matrix reference that will be used to define the preconditioner - * \param[in] geometry_ref - - * \param[in] config_ref - - */ - CJacobiTransposedPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); - - /*! - * \brief destructor of the class - */ - ~CJacobiTransposedPreconditioner() {} - - /*! - * \brief operator that defines the preconditioner operation - * \param[in] u - CSysVector that is being preconditioned - * \param[out] v - CSysVector that is the result of the preconditioning - */ - void operator()(const CSysVector & u, CSysVector & v) const; -}; - -/*! - * \class CILUPreconditioner - * \brief specialization of preconditioner that uses CSysMatrix class - */ -template -class CILUPreconditioner : public CPreconditioner { -private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ - CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ - CConfig* config; /*!< \brief pointer to matrix that defines the config. */ - -public: - - /*! - * \brief constructor of the class - * \param[in] matrix_ref - matrix reference that will be used to define the preconditioner - * \param[in] geometry_ref - - * \param[in] config_ref - - */ - CILUPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); - - /*! - * \brief destructor of the class - */ - ~CILUPreconditioner() {} - - /*! - * \brief operator that defines the preconditioner operation - * \param[in] u - CSysVector that is being preconditioned - * \param[out] v - CSysVector that is the result of the preconditioning - */ - void operator()(const CSysVector & u, CSysVector & v) const; -}; - -/*! - * \class CLU_SGSPreconditioner - * \brief specialization of preconditioner that uses CSysMatrix class - */ -template -class CLU_SGSPreconditioner : public CPreconditioner { -private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ - CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ - CConfig* config; /*!< \brief pointer to matrix that defines the config. */ - -public: - - /*! - * \brief constructor of the class - * \param[in] matrix_ref - matrix reference that will be used to define the preconditioner - * \param[in] geometry_ref - - * \param[in] config_ref - - */ - CLU_SGSPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); - - /*! - * \brief destructor of the class - */ - ~CLU_SGSPreconditioner() {} - - /*! - * \brief operator that defines the preconditioner operation - * \param[in] u - CSysVector that is being preconditioned - * \param[out] v - CSysVector that is the result of the preconditioning - */ - void operator()(const CSysVector & u, CSysVector & v) const; -}; - -/*! - * \class CLineletPreconditioner - * \brief specialization of preconditioner that uses CSysMatrix class - */ -template -class CLineletPreconditioner : public CPreconditioner { -private: - CSysMatrix* sparse_matrix; /*!< \brief pointer to matrix that defines the preconditioner. */ - CGeometry* geometry; /*!< \brief pointer to matrix that defines the geometry. */ - CConfig* config; /*!< \brief pointer to matrix that defines the config. */ - -public: - - /*! - * \brief constructor of the class - * \param[in] matrix_ref - matrix reference that will be used to define the preconditioner - * \param[in] geometry_ref - - * \param[in] config_ref - - */ - CLineletPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref); - - /*! - * \brief destructor of the class - */ - ~CLineletPreconditioner() {} - - /*! - * \brief operator that defines the preconditioner operation - * \param[in] u - CSysVector that is being preconditioned - * \param[out] v - CSysVector that is the result of the preconditioning - */ - void operator()(const CSysVector & u, CSysVector & v) const; -}; - -#include "matrix_structure.inl" +template<> template<> +inline passivedouble CSysMatrix::ActiveAssign(const su2double & val) const { return SU2_TYPE::GetValue(val); } +#endif diff --git a/Common/include/matrix_structure.inl b/Common/include/matrix_structure.inl index 1d02bad80b6d..e800fbaf31ae 100644 --- a/Common/include/matrix_structure.inl +++ b/Common/include/matrix_structure.inl @@ -1,6 +1,10 @@ /*! * \file matrix_structure.inl * \brief In-Line subroutines of the matrix_structure.hpp file. + * \note These are the "private" inlines, they are not needed outside of + * the .cpp file and so they are hidden to avoid triggering recompilation + * of other units when changes are made here. + * * \author F. Palacios, A. Bueno, T. Economon * \version 6.2.0 "Falcon" * @@ -37,48 +41,27 @@ #pragma once -template -inline void CSysMatrix::SetValZero(void) { - if(NULL != matrix) { - for (unsigned long index = 0; index < nnz*nVar*nEqn; index++) - matrix[index] = 0.0; - } -} +#include "matrix_structure.hpp" template -template -inline DstType CSysMatrix::ActiveAssign(const SrcType & val) const { return val; } +inline ScalarType *CSysMatrix::GetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j) { -#ifdef CODI_REVERSE_TYPE -template<> template<> -inline passivedouble CSysMatrix::ActiveAssign(const su2double & val) const { return SU2_TYPE::GetValue(val); } - -template<> template<> -inline passivedouble CSysMatrix::ActiveAssign(const su2double & val) const { return SU2_TYPE::GetValue(val); } -#endif + for (unsigned long index = row_ptr_ilu[block_i]; index < row_ptr_ilu[block_i+1]; index++) + if (col_ind_ilu[index] == block_j) + return &(ILU_matrix[index*nVar*nEqn]); -template -template -inline DstType CSysMatrix::PassiveAssign(const SrcType & val) const { -#if defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE) - return SU2_TYPE::GetValue(val); -#else - return val; -#endif + return NULL; } template -template -inline void CSysMatrix::SetBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block) { +inline void CSysMatrix::SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block) { - unsigned long iVar, jVar, index, step = 0; + unsigned long iVar, index; - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) - matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = PassiveAssign(val_block[iVar][jVar]); + for (index = row_ptr_ilu[block_i]; index < row_ptr_ilu[block_i+1]; index++) { + if (col_ind_ilu[index] == block_j) { + for (iVar = 0; iVar < nVar*nEqn; iVar++) + ILU_matrix[index*nVar*nEqn+iVar] = val_block[iVar]; break; } } @@ -86,197 +69,291 @@ inline void CSysMatrix::SetBlock(unsigned long block_i, unsigned lon } template -template -inline void CSysMatrix::SetBlock(unsigned long block_i, unsigned long block_j, OtherType *val_block) { - - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { +inline void CSysMatrix::SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block) { + + unsigned long iVar, jVar, index; + + for (index = row_ptr_ilu[block_i]; index < row_ptr_ilu[block_i+1]; index++) { + if (col_ind_ilu[index] == block_j) { for (iVar = 0; iVar < nVar; iVar++) for (jVar = 0; jVar < nEqn; jVar++) - matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = PassiveAssign(val_block[iVar*nVar+jVar]); + ILU_matrix[index*nVar*nEqn+iVar*nEqn+jVar] = val_block[jVar*nVar+iVar]; break; } } - + } template -template -inline void CSysMatrix::AddBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block) { +inline void CSysMatrix::SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block) { - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) - matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] += PassiveAssign(val_block[iVar][jVar]); + for (unsigned long index = row_ptr_ilu[block_i]; index < row_ptr_ilu[block_i+1]; index++) { + if (col_ind_ilu[index] == block_j) { + MatrixSubtraction(&ILU_matrix[index*nVar*nEqn], val_block, &ILU_matrix[index*nVar*nEqn]); break; } } } -template -template -inline void CSysMatrix::SubtractBlock(unsigned long block_i, unsigned long block_j, OtherType **val_block) { - - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) - matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] -= PassiveAssign(val_block[iVar][jVar]); - break; +template +inline void gemv_impl(const unsigned long n, const T *a, const T *b, T *c) { + /*--- + This is a templated version of GEMV with the constants as boolean + template parameters so that they can be optimized away at compilation. + This is still the traditional "row dot vector" method. + ---*/ + unsigned long i, j; + for (i = 0; i < n; i++) { + if (!beta) c[i] = 0.0; + for (j = 0; j < n; j++) { + if (alpha) c[transp? j:i] += a[i*n+j] * b[transp? i:j]; + else c[transp? j:i] -= a[i*n+j] * b[transp? i:j]; } } - } -template -template -inline void CSysMatrix::AddVal2Diag(unsigned long block_i, OtherType val_matrix) { - - unsigned long step = 0, iVar, index; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_i) { // Only elements on the diagonal - for (iVar = 0; iVar < nVar; iVar++) - matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+iVar] += PassiveAssign(val_matrix); - break; +template +inline void gemm_impl(const unsigned long n, const T *a, const T *b, T *c) { + /*--- Same deal as for GEMV but here only the type is templated. ---*/ + unsigned long i, j, k; + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + c[i*n+j] = 0.0; + for (k = 0; k < n; k++) + c[i*n+j] += a[i*n+k] * b[k*n+j]; } } - } -template -template -inline void CSysMatrix::SetVal2Diag(unsigned long block_i, OtherType val_matrix) { - - unsigned long step = 0, iVar, jVar, index; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_i) { // Only elements on the diagonal - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+jVar] = 0.0; - - for (iVar = 0; iVar < nVar; iVar++) - matrix[(row_ptr[block_i]+step-1)*nVar*nVar+iVar*nVar+iVar] = PassiveAssign(val_matrix); - - break; - } - } - +#define __MATVECPROD_SIGNATURE__(TYPE,NAME) \ +inline void CSysMatrix::NAME(const TYPE *matrix, const TYPE *vector, TYPE *product) + +#define MATVECPROD_SIGNATURE(NAME) template __MATVECPROD_SIGNATURE__(ScalarType,NAME) + +#if !defined(USE_MKL) +MATVECPROD_SIGNATURE( MatrixVectorProduct ) { + /*--- + Without MKL (default) picture copying the body of gemv_impl + here and resolving the conditionals at compilation. + ---*/ + gemv_impl(nVar, matrix, vector, product); } -template -inline CSysMatrixVectorProduct::CSysMatrixVectorProduct(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { - sparse_matrix = &matrix_ref; - geometry = geometry_ref; - config = config_ref; +MATVECPROD_SIGNATURE( MatrixVectorProductAdd ) { + gemv_impl(nVar, matrix, vector, product); } -template -inline void CSysMatrixVectorProduct::operator()(const CSysVector & u, CSysVector & v) const { - if (sparse_matrix == NULL) { - cerr << "CSysMatrixVectorProduct::operator()(const CSysVector &, CSysVector &): " << endl; - cerr << "pointer to sparse matrix is NULL." << endl; - throw(-1); - } - sparse_matrix->MatrixVectorProduct(u, v, geometry, config); +MATVECPROD_SIGNATURE( MatrixVectorProductSub ) { + gemv_impl(nVar, matrix, vector, product); } -template -inline CSysMatrixVectorProductTransposed::CSysMatrixVectorProductTransposed(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { - sparse_matrix = &matrix_ref; - geometry = geometry_ref; - config = config_ref; +MATVECPROD_SIGNATURE( MatrixVectorProductTransp ) { + gemv_impl(nVar, matrix, vector, product); } template -inline void CSysMatrixVectorProductTransposed::operator()(const CSysVector & u, CSysVector & v) const { - if (sparse_matrix == NULL) { - cerr << "CSysMatrixVectorProduct::operator()(const CSysVector &, CSysVector &): " << endl; - cerr << "pointer to sparse matrix is NULL." << endl; - throw(-1); - } - sparse_matrix->MatrixVectorProductTransposed(u, v, geometry, config); +inline void CSysMatrix::MatrixMatrixProduct(const ScalarType *matrix_a, const ScalarType *matrix_b, ScalarType *product) { + gemm_impl(nVar, matrix_a, matrix_b, product); +} +#else +MATVECPROD_SIGNATURE( MatrixVectorProduct ) { + /*--- With MKL we use the just-in-time kernels instead of the naive implementation. ---*/ + MatrixVectorProductKernelBetaZero(MatrixVectorProductJitterBetaZero, const_cast(vector), + const_cast(matrix), product ); } -template -inline CJacobiPreconditioner::CJacobiPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { - sparse_matrix = &matrix_ref; - geometry = geometry_ref; - config = config_ref; +MATVECPROD_SIGNATURE( MatrixVectorProductAdd ) { + MatrixVectorProductKernelBetaOne(MatrixVectorProductJitterBetaOne, const_cast(vector), + const_cast(matrix), product ); } -template -inline void CJacobiPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { - if (sparse_matrix == NULL) { - cerr << "CJacobiPreconditioner::operator()(const CSysVector &, CSysVector &): " << endl; - cerr << "pointer to sparse matrix is NULL." << endl; - throw(-1); - } - sparse_matrix->ComputeJacobiPreconditioner(u, v, geometry, config); +MATVECPROD_SIGNATURE( MatrixVectorProductSub ) { + MatrixVectorProductKernelAlphaMinusOne(MatrixVectorProductJitterAlphaMinusOne, const_cast(vector), + const_cast(matrix), product ); +} + +MATVECPROD_SIGNATURE( MatrixVectorProductTransp ) { + MatrixVectorProductTranspKernelBetaOne(MatrixVectorProductTranspJitterBetaOne, const_cast(matrix), + const_cast(vector), product ); } template -inline CILUPreconditioner::CILUPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { - sparse_matrix = &matrix_ref; - geometry = geometry_ref; - config = config_ref; +inline void CSysMatrix::MatrixMatrixProduct(const ScalarType *matrix_a, const ScalarType *matrix_b, ScalarType *product) { + MatrixMatrixProductKernel(MatrixMatrixProductJitter, const_cast(matrix_a), + const_cast(matrix_b), product ); +} +#ifdef CODI_REVERSE_TYPE +/*--- WHEN using MKL, AND compiling for AD, we need to specialize for su2double to avoid mixing incompatible types. ---*/ +#define MATVECPROD_SPECIALIZATION(NAME) template<> __MATVECPROD_SIGNATURE__(su2double,NAME) +MATVECPROD_SPECIALIZATION( MatrixVectorProduct ) { + gemv_impl(nVar, matrix, vector, product); +} + +MATVECPROD_SPECIALIZATION( MatrixVectorProductAdd ) { + gemv_impl(nVar, matrix, vector, product); +} + +MATVECPROD_SPECIALIZATION( MatrixVectorProductSub ) { + gemv_impl(nVar, matrix, vector, product); +} + +MATVECPROD_SPECIALIZATION( MatrixVectorProductTransp ) { + gemv_impl(nVar, matrix, vector, product); +} + +template<> +inline void CSysMatrix::MatrixMatrixProduct(const su2double *matrix_a, const su2double *matrix_b, su2double *product) { + gemm_impl(nVar, matrix_a, matrix_b, product); } +#undef MATVECPROD_SPECIALIZATION +#endif // CODI_REVERSE_TYPE +#endif // USE_MKL + +#undef MATVECPROD_SIGNATURE +#undef __MATVECPROD_SIGNATURE__ template -inline void CILUPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { - if (sparse_matrix == NULL) { - cerr << "CILUPreconditioner::operator()(const CSysVector &, CSysVector &): " << endl; - cerr << "pointer to sparse matrix is NULL." << endl; - throw(-1); +inline void CSysMatrix::Gauss_Elimination(ScalarType* matrix, ScalarType* vec) { + + /*--- + This is a relatively large method to inline but maybe better + code will be generated for the special case nVar=1 this way. + ---*/ + + if (nVar==1) {vec[0] /= matrix[0]; return;} + +#ifdef USE_MKL_LAPACK + // With MKL_DIRECT_CALL enabled, this is significantly faster than native code on Intel Architectures. + LAPACKE_dgetrf( LAPACK_ROW_MAJOR, nVar, nVar, matrix, nVar, mkl_ipiv ); + LAPACKE_dgetrs( LAPACK_ROW_MAJOR, 'N', nVar, 1, matrix, nVar, mkl_ipiv, vec, 1 ); +#else + int iVar, jVar, kVar, nvar = int(nVar); + ScalarType weight; + + /*--- Transform system in Upper Matrix ---*/ + for (iVar = 1; iVar < nvar; iVar++) { + for (jVar = 0; jVar < iVar; jVar++) { + weight = matrix[iVar*nvar+jVar] / matrix[jVar*nvar+jVar]; + for (kVar = jVar; kVar < nvar; kVar++) + matrix[iVar*nvar+kVar] -= weight*matrix[jVar*nvar+kVar]; + vec[iVar] -= weight*vec[jVar]; + } } - sparse_matrix->ComputeILUPreconditioner(u, v, geometry, config); + + /*--- Backwards substitution ---*/ + for (iVar = nvar-1; iVar >= 0; iVar--) { + for (jVar = iVar+1; jVar < nvar; jVar++) + vec[iVar] -= matrix[iVar*nvar+jVar]*vec[jVar]; + vec[iVar] /= matrix[iVar*nvar+iVar]; + } +#endif } template -inline CLU_SGSPreconditioner::CLU_SGSPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { - sparse_matrix = &matrix_ref; - geometry = geometry_ref; - config = config_ref; +inline void CSysMatrix::MatrixInverse(const ScalarType *matrix, ScalarType *inverse) { + + /*--- + This is a generalization of Gaussian elimination for multiple rhs' (the basis vectors). + We could call "Gauss_Elimination" multiple times or fully generalize it for multiple rhs, + the performance of both routines would suffer in both cases without the use of exotic templating. + And so it feels reasonable to have some duplication here. + ---*/ + + if (nVar==1) {inverse[0] = 1.0/matrix[0]; return;} + + int iVar, jVar, nvar = int(nVar); + + /*--- Initialize the inverse and make a copy of the matrix ---*/ + for (iVar = 0; iVar < nvar; iVar++) { + for (jVar = 0; jVar < nvar; jVar++) { + block[iVar*nvar+jVar] = matrix[iVar*nvar+jVar]; + inverse[iVar*nvar+jVar] = ScalarType(iVar==jVar); // identity + } + } + + /*--- Inversion ---*/ +#ifdef USE_MKL_LAPACK + // With MKL_DIRECT_CALL enabled, this is significantly faster than native code on Intel Architectures. + LAPACKE_dgetrf( LAPACK_ROW_MAJOR, nVar, nVar, block, nVar, mkl_ipiv ); + LAPACKE_dgetrs( LAPACK_ROW_MAJOR, 'N', nVar, nVar, block, nVar, mkl_ipiv, inverse, nVar ); +#else + int kVar; + ScalarType weight; + + /*--- Transform system in Upper Matrix ---*/ + for (iVar = 1; iVar < nvar; iVar++) { + for (jVar = 0; jVar < iVar; jVar++) + { + weight = block[iVar*nvar+jVar] / block[jVar*nvar+jVar]; + + for (kVar = jVar; kVar < nvar; kVar++) + block[iVar*nvar+kVar] -= weight*block[jVar*nvar+kVar]; + + /*--- at this stage "inverse" is lower triangular so not all cols need updating ---*/ + for (kVar = 0; kVar <= jVar; kVar++) + inverse[iVar*nvar+kVar] -= weight*inverse[jVar*nvar+kVar]; + } + } + + /*--- Backwards substitution ---*/ + for (iVar = nvar-1; iVar >= 0; iVar--) + { + for (jVar = iVar+1; jVar < nvar; jVar++) + for (kVar = 0; kVar < nvar; kVar++) + inverse[iVar*nvar+kVar] -= block[iVar*nvar+jVar] * inverse[jVar*nvar+kVar]; + + for (kVar = 0; kVar < nvar; kVar++) + inverse[iVar*nvar+kVar] /= block[iVar*nvar+iVar]; + } +#endif } template -inline void CLU_SGSPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { - if (sparse_matrix == NULL) { - cerr << "CLU_SGSPreconditioner::operator()(const CSysVector &, CSysVector &): " << endl; - cerr << "pointer to sparse matrix is NULL." << endl; - throw(-1); +inline void CSysMatrix::Gauss_Elimination(unsigned long block_i, ScalarType* rhs, bool transposed) { + + unsigned long iVar, jVar; + ScalarType *Block = GetBlock(block_i, block_i); + + /*--- Copy block, as the algorithm modifies the matrix ---*/ + + if (!transposed) { + // If source and dest overlap higher level problems occur, so memcpy is safe. And it is faster. + memcpy( block, Block, nVar*nVar*sizeof(ScalarType) ); + +// for (iVar = 0; iVar < nVar*nVar; iVar++) +// block[iVar] = Block[iVar]; + + } else { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) + block[iVar*nVar+jVar] = Block[jVar*nVar+iVar]; } - sparse_matrix->ComputeLU_SGSPreconditioner(u, v, geometry, config); + + /*--- Solve system ---*/ + + Gauss_Elimination(block, rhs); + } template -inline CLineletPreconditioner::CLineletPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref) { - sparse_matrix = &matrix_ref; - geometry = geometry_ref; - config = config_ref; +inline void CSysMatrix::InverseDiagonalBlock(unsigned long block_i, ScalarType *invBlock, bool transpose) { + + const ScalarType* mat = GetBlock(block_i, block_i); + MatrixInverse(mat, invBlock); + + if (transpose) // swap off-diag + for (unsigned long iVar = 0; iVar < nVar-1; ++iVar) + for (unsigned long jVar = iVar+1; jVar < nVar; ++jVar) { + ScalarType tmp = invBlock[iVar*nVar+jVar]; + invBlock[iVar*nVar+jVar] = invBlock[jVar*nVar+iVar]; + invBlock[jVar*nVar+iVar] = tmp; + } } template -inline void CLineletPreconditioner::operator()(const CSysVector & u, CSysVector & v) const { - if (sparse_matrix == NULL) { - cerr << "CLineletPreconditioner::operator()(const CSysVector &, CSysVector &): " << endl; - cerr << "pointer to sparse matrix is NULL." << endl; - throw(-1); - } - sparse_matrix->ComputeLineletPreconditioner(u, v, geometry, config); +inline void CSysMatrix::InverseDiagonalBlock_ILUMatrix(unsigned long block_i, ScalarType *invBlock) { + + const ScalarType* mat = GetBlock_ILUMatrix(block_i, block_i); + MatrixInverse(mat, invBlock); } diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 6d422ee2e2b3..0037469be6e5 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1700,10 +1700,7 @@ enum ENUM_LINEAR_SOLVER { FGMRES = 5, /*!< \brief Flexible Generalized Minimal Residual method. */ BCGSTAB = 6, /*!< \brief BCGSTAB - Biconjugate Gradient Stabilized Method (main solver). */ RESTARTED_FGMRES = 7, /*!< \brief Flexible Generalized Minimal Residual method with restart. */ - SMOOTHER_LUSGS = 8, /*!< \brief LU_SGS smoother. */ - SMOOTHER_JACOBI = 9, /*!< \brief Jacobi smoother. */ - SMOOTHER_ILU = 10, /*!< \brief ILU smoother. */ - SMOOTHER_LINELET = 11 /*!< \brief Linelet smoother. */ + SMOOTHER = 8, /*!< \brief Iterative smoother. */ }; static const map Linear_Solver_Map = CCreateMap ("STEEPEST_DESCENT", STEEPEST_DESCENT) @@ -1713,10 +1710,7 @@ static const map Linear_Solver_Map = CCreateMap CSysVector operator*(const ScalarType & val, const CSysVector & u); - -/*! - * \class CMatrixVectorProduct - * \brief abstract base class for defining matrix-vector products - * \author J. Hicken. - * - * The Krylov-subspace solvers require only matrix-vector products and - * not the actual matrix/Jacobian. We need some way to indicate which - * function will perform the product. However, sometimes the - * functions that define the product will require different numbers - * and types of inputs. For example, the forward-difference - * approximation to a Jacobian-vector product requires the vector that - * defines the Jacobian and a perturbation parameter. The - * CMatrixVectorProduct class is used to derive child classes that can - * handle the different types of matrix-vector products and still be - * passed to a single implementation of the Krylov solvers. - */ -template -class CMatrixVectorProduct { -public: - virtual ~CMatrixVectorProduct() = 0; ///< class destructor - virtual void operator()(const CSysVector & u, CSysVector & v) - const = 0; ///< matrix-vector product operation -}; -template -inline CMatrixVectorProduct::~CMatrixVectorProduct() {} - -/*! - * \class CPreconditioner - * \brief abstract base class for defining preconditioning operation - * \author J. Hicken. - * - * See the remarks regarding the CMatrixVectorProduct class. The same - * idea applies here to the preconditioning operation. - */ -template -class CPreconditioner { -public: - virtual ~CPreconditioner() = 0; ///< class destructor - virtual void operator()(const CSysVector & u, CSysVector & v) - const = 0; ///< preconditioning operation -}; -template -inline CPreconditioner::~CPreconditioner() {} - -#include "vector_structure.inl" diff --git a/Common/include/vector_structure.inl b/Common/include/vector_structure.inl deleted file mode 100644 index 6d724333ee16..000000000000 --- a/Common/include/vector_structure.inl +++ /dev/null @@ -1,108 +0,0 @@ -/*! - * \file vector_structure.inl - * \brief inline subroutines of the vector_structure.hpp file. - * \author F. Palacios, J. Hicken - * \version 6.2.0 "Falcon" - * - * The current SU2 release has been coordinated by the - * SU2 International Developers Society - * with selected contributions from the open-source community. - * - * The main research teams contributing to the current release are: - * - Prof. Juan J. Alonso's group at Stanford University. - * - Prof. Piero Colonna's group at Delft University of Technology. - * - Prof. Nicolas R. Gauger's group at Kaiserslautern University of Technology. - * - Prof. Alberto Guardone's group at Polytechnic University of Milan. - * - Prof. Rafael Palacios' group at Imperial College London. - * - Prof. Vincent Terrapon's group at the University of Liege. - * - Prof. Edwin van der Weide's group at the University of Twente. - * - Lab. of New Concepts in Aeronautics at Tech. Institute of Aeronautics. - * - * Copyright 2012-2019, Francisco D. Palacios, Thomas D. Economon, - * Tim Albring, and the SU2 contributors. - * - * 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 . - */ - -#pragma once - -template -inline void CSysVector::SetValZero(void) { - for (unsigned long i = 0; i < nElm; i++) - vec_val[i] = 0.0; -} - -template -inline unsigned long CSysVector::GetLocSize() const { return nElm; } - -template -unsigned long CSysVector::GetNElmDomain() const { return nElmDomain; } - -template -inline unsigned long CSysVector::GetSize() const { -#ifdef HAVE_MPI - return nElmGlobal; -#else - return (unsigned long)nElm; -#endif -} - -template -inline unsigned short CSysVector::GetNVar() const { return nVar; } - -template -inline unsigned long CSysVector::GetNBlk() const { return nBlk; } - -template -inline unsigned long CSysVector::GetNBlkDomain() const { return nBlkDomain; } - -template -inline ScalarType & CSysVector::operator[](const unsigned long & i) { return vec_val[i]; } - -template -inline const ScalarType & CSysVector::operator[](const unsigned long & i) const { return vec_val[i]; } - -template -template -void CSysVector::PassiveCopy(const CSysVector& other) { - - /*--- This is a method and not the overload of an operator to make sure who - calls it knows the consequence to the derivative information (lost) ---*/ - - /*--- check if self-assignment, otherwise perform deep copy ---*/ - if ((const void*)this == (const void*)&other) return; - - /*--- determine if (re-)allocation is needed ---*/ - if (nElm != other.GetLocSize() && vec_val != NULL) { - delete [] vec_val; - vec_val = NULL; - } - - /*--- copy ---*/ - nElm = other.GetLocSize(); - nElmDomain = other.GetNElmDomain(); - nBlk = other.GetNBlk(); - nBlkDomain = other.GetNBlkDomain(); - nVar = other.GetNVar(); - - if (vec_val == NULL) - vec_val = new ScalarType[nElm]; - - for (unsigned long i = 0; i < nElm; i++) - vec_val[i] = SU2_TYPE::GetValue(other[i]); - -#ifdef HAVE_MPI - nElmGlobal = other.GetSize(); -#endif -} diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp index c471d6a8d393..3da5f7e6555e 100644 --- a/Common/src/config_structure.cpp +++ b/Common/src/config_structure.cpp @@ -1296,6 +1296,8 @@ void CConfig::SetConfig_Options(unsigned short val_iZone, unsigned short val_nZo addUnsignedShortOption("LINEAR_SOLVER_ILU_FILL_IN", Linear_Solver_ILU_n, 0); /* DESCRIPTION: Maximum number of iterations of the linear solver for the implicit formulation */ addUnsignedLongOption("LINEAR_SOLVER_RESTART_FREQUENCY", Linear_Solver_Restart_Frequency, 10); + /* DESCRIPTION: Relaxation factor for iterative linear smoothers (SMOOTHER_ILU/JACOBI/LU-SGS/LINELET) */ + addDoubleOption("LINEAR_SOLVER_SMOOTHER_RELAXATION", Linear_Solver_Smoother_Relaxation, 1.0); /* DESCRIPTION: Relaxation of the flow equations solver for the implicit formulation */ addDoubleOption("RELAXATION_FACTOR_FLOW", Relaxation_Factor_Flow, 1.0); /* DESCRIPTION: Relaxation of the turb equations solver for the implicit formulation */ @@ -5862,41 +5864,31 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) { cout << "Euler implicit method for the flow equations." << endl; switch (Kind_Linear_Solver) { case BCGSTAB: - cout << "BCGSTAB is used for solving the linear system." << endl; - switch (Kind_Linear_Solver_Prec) { - case ILU: cout << "Using a ILU("<< Linear_Solver_ILU_n <<") preconditioning."<< endl; break; - case LINELET: cout << "Using a linelet preconditioning."<< endl; break; - case LU_SGS: cout << "Using a LU-SGS preconditioning."<< endl; break; - case JACOBI: cout << "Using a Jacobi preconditioning."<< endl; break; - } - cout << "Convergence criteria of the linear solver: "<< Linear_Solver_Error <<"."<< endl; - cout << "Max number of linear iterations: "<< Linear_Solver_Iter <<"."<< endl; - break; case FGMRES: case RESTARTED_FGMRES: - cout << "FGMRES is used for solving the linear system." << endl; + if (Kind_Linear_Solver == BCGSTAB) + cout << "BCGSTAB is used for solving the linear system." << endl; + else + cout << "FGMRES is used for solving the linear system." << endl; switch (Kind_Linear_Solver_Prec) { case ILU: cout << "Using a ILU("<< Linear_Solver_ILU_n <<") preconditioning."<< endl; break; case LINELET: cout << "Using a linelet preconditioning."<< endl; break; - case LU_SGS: cout << "Using a LU-SGS preconditioning."<< endl; break; - case JACOBI: cout << "Using a Jacobi preconditioning."<< endl; break; + case LU_SGS: cout << "Using a LU-SGS preconditioning."<< endl; break; + case JACOBI: cout << "Using a Jacobi preconditioning."<< endl; break; } - cout << "Convergence criteria of the linear solver: "<< Linear_Solver_Error <<"."<< endl; - cout << "Max number of linear iterations: "<< Linear_Solver_Iter <<"."<< endl; - break; - case SMOOTHER_JACOBI: - cout << "A Jacobi method is used for smoothing the linear system." << endl; - break; - case SMOOTHER_ILU: - cout << "A ILU("<< Linear_Solver_ILU_n <<") method is used for smoothing the linear system." << endl; break; - case SMOOTHER_LUSGS: - cout << "A LU-SGS method is used for smoothing the linear system." << endl; - break; - case SMOOTHER_LINELET: - cout << "A Linelet method is used for smoothing the linear system." << endl; + case SMOOTHER: + switch (Kind_Linear_Solver_Prec) { + case ILU: cout << "A ILU(" << Linear_Solver_ILU_n << ")"; break; + case LINELET: cout << "A Linelet"; break; + case LU_SGS: cout << "A LU-SGS"; break; + case JACOBI: cout << "A Jacobi"; break; + } + cout << " method is used for smoothing the linear system." << endl; break; } + cout << "Convergence criteria of the linear solver: "<< Linear_Solver_Error <<"."<< endl; + cout << "Max number of linear iterations: "<< Linear_Solver_Iter <<"."<< endl; break; case CLASSICAL_RK4_EXPLICIT: cout << "Classical RK4 explicit method for the flow equations." << endl; diff --git a/Common/src/grid_movement_structure.cpp b/Common/src/grid_movement_structure.cpp index eefb3d0fc733..2734055c89c4 100644 --- a/Common/src/grid_movement_structure.cpp +++ b/Common/src/grid_movement_structure.cpp @@ -9528,130 +9528,35 @@ void CElasticityMovement::SetBoundaryDisplacements(CGeometry *geometry, CConfig void CElasticityMovement::SetClamped_Boundary(CGeometry *geometry, CConfig *config, unsigned short val_marker){ unsigned long iNode, iVertex; - unsigned long iPoint, jPoint; + + Solution[0] = 0.0; + Solution[1] = 0.0; + if (nDim==3) Solution[2] = 0.0; for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Get node index ---*/ - iNode = geometry->vertex[val_marker][iVertex]->GetNode(); - if (geometry->node[iNode]->GetDomain()) { - - if (nDim == 2) { - Solution[0] = 0.0; Solution[1] = 0.0; - Residual[0] = 0.0; Residual[1] = 0.0; - } - else { - Solution[0] = 0.0; Solution[1] = 0.0; Solution[2] = 0.0; - Residual[0] = 0.0; Residual[1] = 0.0; Residual[2] = 0.0; - } - - /*--- Initialize the reaction vector ---*/ - - LinSysRes.SetBlock(iNode, Residual); - LinSysSol.SetBlock(iNode, Solution); - - /*--- STRONG ENFORCEMENT OF THE CLAMPED BOUNDARY CONDITION ---*/ - - /*--- Delete the full row for node iNode ---*/ - for (jPoint = 0; jPoint < nPoint; jPoint++){ - if (iNode != jPoint) { - StiffMatrix.SetBlock(iNode,jPoint,matrixZeros); - } - else{ - StiffMatrix.SetBlock(iNode,jPoint,matrixId); - } - } + /*--- Set and enforce solution ---*/ + LinSysSol.SetBlock(iNode, Solution); + StiffMatrix.EnforceSolutionAtNode(iNode, Solution, LinSysRes); - /*--- Delete the full column for node iNode ---*/ - for (iPoint = 0; iPoint < nPoint; iPoint++){ - if (iNode != iPoint) { - StiffMatrix.SetBlock(iPoint,iNode,matrixZeros); - } - } - } - else { - /*--- Delete the column (iPoint is halo so Send/Recv does the rest) ---*/ - for (iPoint = 0; iPoint < nPoint; iPoint++) StiffMatrix.SetBlock(iPoint,iNode,matrixZeros); - } } } void CElasticityMovement::SetMoving_Boundary(CGeometry *geometry, CConfig *config, unsigned short val_marker){ - unsigned short iDim, jDim; - unsigned long iNode, iVertex; - unsigned long iPoint, jPoint; - - su2double valJacobian_ij_00 = 0.0; - su2double auxJacobian_ij[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Get node index ---*/ - iNode = geometry->vertex[val_marker][iVertex]->GetNode(); - /*--- Get the displacement of the node ---*/ - - for (iDim = 0; iDim < nDim; iDim++) - Solution[iDim] = LinSysSol.GetBlock(iNode, iDim); - - if (geometry->node[iNode]->GetDomain()) { - - /*--- Initialize the reaction vector ---*/ - - LinSysRes.SetBlock(iNode, Solution); - - /*--- STRONG ENFORCEMENT OF THE DISPLACEMENT BOUNDARY CONDITION ---*/ - - /*--- Delete the full row for node iNode ---*/ - for (jPoint = 0; jPoint < nPoint; jPoint++){ - if (iNode != jPoint) { - StiffMatrix.SetBlock(iNode,jPoint,matrixZeros); - } - else{ - StiffMatrix.SetBlock(iNode,jPoint,matrixId); - } - } - } - - /*--- Always delete the iNode column, even for halos ---*/ - for (iPoint = 0; iPoint < nPoint; iPoint++){ - - /*--- Check if the term K(iPoint, iNode) is 0 ---*/ - valJacobian_ij_00 = StiffMatrix.GetBlock(iPoint,iNode,0,0); - - /*--- If the node iNode has a crossed dependency with the point iPoint ---*/ - if (valJacobian_ij_00 != 0.0 ){ - - /*--- Retrieve the Jacobian term ---*/ - for (iDim = 0; iDim < nDim; iDim++){ - for (jDim = 0; jDim < nDim; jDim++){ - auxJacobian_ij[iDim][jDim] = StiffMatrix.GetBlock(iPoint,iNode,iDim,jDim); - } - } - - /*--- Multiply by the imposed displacement ---*/ - for (iDim = 0; iDim < nDim; iDim++){ - Residual[iDim] = 0.0; - for (jDim = 0; jDim < nDim; jDim++){ - Residual[iDim] += auxJacobian_ij[iDim][jDim] * Solution[jDim]; - } - } - - /*--- For the whole column, except the diagonal term ---*/ - if (iNode != iPoint) { - /*--- The term is substracted from the residual (right hand side) ---*/ - LinSysRes.SubtractBlock(iPoint, Residual); - /*--- The Jacobian term is now set to 0 ---*/ - StiffMatrix.SetBlock(iPoint,iNode,matrixZeros); - } - } - } + /*--- Enforce solution ---*/ + StiffMatrix.EnforceSolutionAtNode(iNode, LinSysSol.GetBlock(iNode), LinSysRes); } diff --git a/Common/src/interpolation_structure.cpp b/Common/src/interpolation_structure.cpp index 3f0d1a81ccc0..a0798673959b 100644 --- a/Common/src/interpolation_structure.cpp +++ b/Common/src/interpolation_structure.cpp @@ -37,6 +37,19 @@ #include "../include/interpolation_structure.hpp" +#if defined(HAVE_MKL) +#include "mkl.h" +#ifndef HAVE_LAPACK +#define HAVE_LAPACK +#endif +#elif defined(HAVE_LAPACK) +/*--- Lapack / Blas routines used in RBF interpolation. ---*/ +extern "C" void dsptrf_(char*, int*, passivedouble*, int*, int*); +extern "C" void dsptri_(char*, int*, passivedouble*, int*, passivedouble*, int*); +extern "C" void dsymm_(char*, char*, int*, int*, passivedouble*, passivedouble*, int*, + passivedouble*, int*, passivedouble*, passivedouble*, int*); +#endif + CInterpolator::CInterpolator(void) { size = SU2_MPI::GetSize(); @@ -3687,7 +3700,7 @@ void CSymmetricMatrix::Invert(const bool is_spd) if(!is_spd) LUDecompose(); else CholeskyDecompose(true); CalcInv(true); -#endif // HAVE_LAPACK +#endif } void CSymmetricMatrix::MatVecMult(passivedouble *v) @@ -3722,9 +3735,7 @@ void CSymmetricMatrix::MatMatMult(bool left_side, su2double *mat_vec_in, int N) } #ifdef HAVE_LAPACK - - char side[1]={'R'}, uplo[1]={'L'}; // Right side because mat_vec in row major order - passivedouble *val_full, alpha=1, beta=0; + passivedouble *val_full, alpha=1.0, beta=0.0; /*--- Copy packed storage to full storage to use BLAS level 3 routine ---*/ val_full = new passivedouble [sz*sz]; @@ -3737,7 +3748,12 @@ void CSymmetricMatrix::MatMatMult(bool left_side, su2double *mat_vec_in, int N) } } +#ifndef HAVE_MKL + char side[1]={'R'}, uplo[1]={'L'}; // Right side because mat_vec in row major order dsymm_(side, uplo, &N, &sz, &alpha, val_full, &sz, mat_vec, &N, &beta, tmp_res, &N); +#else + cblas_dsymm(CblasRowMajor, CblasRight, CblasLower, N, sz, alpha, val_full, sz, mat_vec, N, beta, tmp_res, N); +#endif delete [] val_full; diff --git a/Common/src/linear_solvers_structure.cpp b/Common/src/linear_solvers_structure.cpp index 166cfec38c9e..9409f0829145 100644 --- a/Common/src/linear_solvers_structure.cpp +++ b/Common/src/linear_solvers_structure.cpp @@ -39,8 +39,8 @@ #include "../include/linear_solvers_structure_b.hpp" template -CSysSolve::CSysSolve(const bool mesh_deform_mode) : cg_ready(false), bcg_ready(false), gmres_ready(false) { - +CSysSolve::CSysSolve(const bool mesh_deform_mode) : cg_ready(false), bcg_ready(false), + gmres_ready(false), smooth_ready(false) { mesh_deform = mesh_deform_mode; LinSysRes_ptr = NULL; LinSysSol_ptr = NULL; @@ -396,6 +396,7 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector::BCGSTAB_LinSolver(const CSysVector +unsigned long CSysSolve::Smoother_LinSolver(const CSysVector & b, CSysVector & x, + CMatrixVectorProduct & mat_vec, CPreconditioner & precond, + ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CConfig *config) { + + int rank = SU2_MPI::GetRank(); + ScalarType norm_r = 0.0, norm0 = 0.0; + unsigned long i = 0; + + /*--- Relaxation factor, see comments inside the loop over the smoothing iterations. ---*/ + ScalarType omega = SU2_TYPE::GetValue(config->GetLinear_Solver_Smoother_Relaxation()); + + if (m < 1) { + char buf[100]; + SPRINTF(buf, "Illegal value for smoothing iterations, m = %lu", m ); + SU2_MPI::Error(string(buf), CURRENT_FUNCTION); + } + + /*--- Allocate vectors for residual (r), solution increment (z), and matrix-vector + product (A_x), for the latter two this is done only on the first call to the method. ---*/ + + if (!smooth_ready) { + z = b; + A_x = b; + smooth_ready = true; + } + + /*--- Compute the initial residual and check if the system is already solved (if in COMM_FULL mode). ---*/ + + mat_vec(x, A_x); + r = b; r -= A_x; + + /*--- Only compute the residuals in full communication mode. ---*/ + + if (config->GetComm_Level() == COMM_FULL) { + + norm_r = r.norm(); + norm0 = b.norm(); + if ( (norm_r < tol*norm0) || (norm_r < eps) ) { + if (rank == MASTER_NODE) cout << "CSysSolve::Smoother_LinSolver(): system solved by initial guess." << endl; + return 0; + } + + /*--- Set the norm to the initial initial residual value. ---*/ + + norm0 = norm_r; + + /*--- Output header information including initial residual. ---*/ + + if ((monitoring) && (rank == MASTER_NODE)) { + cout << "\n# " << "Smoother" << " residual history" << endl; + cout << "# Residual tolerance target = " << tol << endl; + cout << "# Initial residual norm = " << norm_r << endl; + cout << " " << i << " " << norm_r/norm0 << endl; + } + + } + + /*--- Smoothing Iterations ---*/ + + for (i=0; i1 is NOT equivalent to SOR. ---*/ + + x.Plus_AX(omega, z); + r.Plus_AX(-omega, A_x); + + /*--- Only compute the residuals in full communication mode. ---*/ + /*--- Check if solution has converged, else output the relative residual if necessary. ---*/ + + if (config->GetComm_Level() == COMM_FULL) { + norm_r = r.norm(); + if (norm_r < tol*norm0) break; + if (((monitoring) && (rank == MASTER_NODE)) && ((i+1) % 5 == 0)) + cout << " " << i << " " << norm_r/norm0 << endl; + } + } + + if ((monitoring) && (rank == MASTER_NODE) && (config->GetComm_Level() == COMM_FULL)) { + cout << "# Smoother final (true) residual:" << endl; + cout << "# Iteration = " << i << ": |res|/|res0| = " << norm_r/norm0 << ".\n" << endl; + } + + (*residual) = norm_r; + return i; +} + template<> void CSysSolve::HandleTemporariesIn(CSysVector & LinSysRes, CSysVector & LinSysSol) { @@ -738,7 +842,7 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CS ScreenOutput = config->GetDeform_Output(); } - CMatrixVectorProduct *mat_vec = NULL; + /*--- Stop the recording for the linear solver ---*/ bool TapeActive = NO; @@ -751,102 +855,69 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CS AD::SetExtFuncIn(&LinSysRes[0], LinSysRes.GetLocSize()); - /*--- Stop the recording for the linear solver ---*/ - AD::StopRecording(); #endif } - /*--- Solve the linear system using a Krylov subspace method ---*/ + /*--- Create matrix-vector product, preconditioner, and solve the linear system ---*/ HandleTemporariesIn(LinSysRes, LinSysSol); - if (KindSolver == BCGSTAB || KindSolver == CONJUGATE_GRADIENT || - KindSolver == FGMRES || KindSolver == RESTARTED_FGMRES ) { - - mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); - CPreconditioner* precond = NULL; - - switch (KindPrecond) { - case JACOBI: - Jacobian.BuildJacobiPreconditioner(); - precond = new CJacobiPreconditioner(Jacobian, geometry, config); - break; - case ILU: - Jacobian.BuildILUPreconditioner(); - precond = new CILUPreconditioner(Jacobian, geometry, config); - break; - case LU_SGS: - precond = new CLU_SGSPreconditioner(Jacobian, geometry, config); - break; - case LINELET: - Jacobian.BuildJacobiPreconditioner(); - precond = new CLineletPreconditioner(Jacobian, geometry, config); - break; - default: - Jacobian.BuildJacobiPreconditioner(); - precond = new CJacobiPreconditioner(Jacobian, geometry, config); - break; - } - - switch (KindSolver) { - case BCGSTAB: - IterLinSol = BCGSTAB_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); - break; - case FGMRES: - IterLinSol = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); - break; - case CONJUGATE_GRADIENT: - IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); - break; - case RESTARTED_FGMRES: - IterLinSol = 0; - Norm0 = LinSysRes_ptr->norm(); - while (IterLinSol < MaxIter) { - /*--- Enforce a hard limit on total number of iterations ---*/ - unsigned long IterLimit = min(RestartIter, MaxIter-IterLinSol); - IterLinSol += FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, IterLimit, &Residual, ScreenOutput, config); - if ( Residual < SolverTol*Norm0 ) break; - } - break; - } - - /*--- Dealocate memory of the Krylov subspace method ---*/ - - delete mat_vec; - delete precond; - + CMatrixVectorProduct* mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); + CPreconditioner* precond = NULL; + + switch (KindPrecond) { + case JACOBI: + Jacobian.BuildJacobiPreconditioner(); + precond = new CJacobiPreconditioner(Jacobian, geometry, config); + break; + case ILU: + Jacobian.BuildILUPreconditioner(); + precond = new CILUPreconditioner(Jacobian, geometry, config); + break; + case LU_SGS: + precond = new CLU_SGSPreconditioner(Jacobian, geometry, config); + break; + case LINELET: + Jacobian.BuildJacobiPreconditioner(); + precond = new CLineletPreconditioner(Jacobian, geometry, config); + break; + default: + Jacobian.BuildJacobiPreconditioner(); + precond = new CJacobiPreconditioner(Jacobian, geometry, config); + break; } - /*--- Smooth the linear system. ---*/ - - else { - switch (KindSolver) { - case SMOOTHER_LUSGS: - mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); - IterLinSol = Jacobian.LU_SGS_Smoother(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, SolverTol, MaxIter, &Residual, ScreenOutput, geometry, config); - delete mat_vec; - break; - case SMOOTHER_JACOBI: - mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); - Jacobian.BuildJacobiPreconditioner(); - IterLinSol = Jacobian.Jacobi_Smoother(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, SolverTol, MaxIter, &Residual, ScreenOutput, geometry, config); - delete mat_vec; - break; - case SMOOTHER_ILU: - mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config); - Jacobian.BuildILUPreconditioner(); - IterLinSol = Jacobian.ILU_Smoother(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, SolverTol, MaxIter, &Residual, ScreenOutput, geometry, config); - delete mat_vec; - break; - case SMOOTHER_LINELET: - Jacobian.BuildJacobiPreconditioner(); - Jacobian.ComputeLineletPreconditioner(*LinSysRes_ptr, *LinSysSol_ptr, geometry, config); - IterLinSol = 1; - break; - } + switch (KindSolver) { + case BCGSTAB: + IterLinSol = BCGSTAB_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); + break; + case FGMRES: + IterLinSol = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); + break; + case CONJUGATE_GRADIENT: + IterLinSol = CG_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); + break; + case RESTARTED_FGMRES: + IterLinSol = 0; + Norm0 = LinSysRes_ptr->norm(); + while (IterLinSol < MaxIter) { + /*--- Enforce a hard limit on total number of iterations ---*/ + unsigned long IterLimit = min(RestartIter, MaxIter-IterLinSol); + IterLinSol += FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, IterLimit, &Residual, ScreenOutput, config); + if ( Residual < SolverTol*Norm0 ) break; + } + break; + case SMOOTHER: + IterLinSol = Smoother_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, *mat_vec, *precond, SolverTol, MaxIter, &Residual, ScreenOutput, config); + break; + default: + SU2_MPI::Error("Unknown type of linear solver.",CURRENT_FUNCTION); } + delete mat_vec; + delete precond; + HandleTemporariesOut(LinSysSol); if(TapeActive) { diff --git a/Common/src/matrix_structure.cpp b/Common/src/matrix_structure.cpp index 4c0e4312ffba..a6cc41d663da 100644 --- a/Common/src/matrix_structure.cpp +++ b/Common/src/matrix_structure.cpp @@ -35,7 +35,7 @@ * License along with SU2. If not, see . */ -#include "../include/matrix_structure.hpp" +#include "../include/matrix_structure.inl" template CSysMatrix::CSysMatrix(void) { @@ -54,7 +54,6 @@ CSysMatrix::CSysMatrix(void) { row_ptr_ilu = NULL; col_ind_ilu = NULL; block = NULL; - prod_block_vector = NULL; prod_row_vector = NULL; aux_vector = NULL; sum_vector = NULL; @@ -62,34 +61,19 @@ CSysMatrix::CSysMatrix(void) { block_weight = NULL; block_inverse = NULL; - /*--- Linelet preconditioner ---*/ - - LineletBool = NULL; - LineletPoint = NULL; - UBlock = NULL; - invUBlock = NULL; - LBlock = NULL; - yVector = NULL; - zVector = NULL; - rVector = NULL; - LFBlock = NULL; - LyVector = NULL; - FzVector = NULL; - max_nElem = 0; - -#if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) - MatrixMatrixProductJitter = NULL; - MatrixVectorProductJitterBetaOne = NULL; - MatrixVectorProductJitterBetaZero = NULL; - useMKL = false; +#ifdef USE_MKL + MatrixMatrixProductJitter = NULL; + MatrixVectorProductJitterBetaOne = NULL; + MatrixVectorProductJitterBetaZero = NULL; + MatrixVectorProductJitterAlphaMinusOne = NULL; + MatrixVectorProductTranspJitterBetaOne = NULL; + mkl_ipiv = NULL; #endif } template CSysMatrix::~CSysMatrix(void) { - - unsigned long iElem; /*--- Memory deallocation ---*/ @@ -107,37 +91,18 @@ CSysMatrix::~CSysMatrix(void) { if (block_weight != NULL) delete [] block_weight; if (block_inverse != NULL) delete [] block_inverse; - if (prod_block_vector != NULL) delete [] prod_block_vector; if (prod_row_vector != NULL) delete [] prod_row_vector; if (aux_vector != NULL) delete [] aux_vector; if (sum_vector != NULL) delete [] sum_vector; if (invM != NULL) delete [] invM; - if (LineletBool != NULL) delete [] LineletBool; - if (LineletPoint != NULL) delete [] LineletPoint; - - for (iElem = 0; iElem < max_nElem; iElem++) { - if (UBlock[iElem] != NULL) delete [] UBlock[iElem]; - if (invUBlock[iElem] != NULL) delete [] invUBlock[iElem]; - if (LBlock[iElem] != NULL) delete [] LBlock[iElem]; - if (yVector[iElem] != NULL) delete [] yVector[iElem]; - if (zVector[iElem] != NULL) delete [] zVector[iElem]; - if (rVector[iElem] != NULL) delete [] rVector[iElem]; - } - if (UBlock != NULL) delete [] UBlock; - if (invUBlock != NULL) delete [] invUBlock; - if (LBlock != NULL) delete [] LBlock; - if (yVector != NULL) delete [] yVector; - if (zVector != NULL) delete [] zVector; - if (rVector != NULL) delete [] rVector; - - if (LFBlock != NULL) delete [] LFBlock; - if (LyVector != NULL) delete [] LyVector; - if (FzVector != NULL) delete [] FzVector; - -#if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) - if ( MatrixMatrixProductJitter != NULL ) mkl_jit_destroy( MatrixMatrixProductJitter ); - if ( MatrixVectorProductJitterBetaZero != NULL ) mkl_jit_destroy( MatrixVectorProductJitterBetaZero ); - if ( MatrixVectorProductJitterBetaOne != NULL ) mkl_jit_destroy( MatrixVectorProductJitterBetaOne ); + +#ifdef USE_MKL + if ( MatrixMatrixProductJitter != NULL ) mkl_jit_destroy( MatrixMatrixProductJitter ); + if ( MatrixVectorProductJitterBetaZero != NULL ) mkl_jit_destroy( MatrixVectorProductJitterBetaZero ); + if ( MatrixVectorProductJitterBetaOne != NULL ) mkl_jit_destroy( MatrixVectorProductJitterBetaOne ); + if ( MatrixVectorProductJitterAlphaMinusOne != NULL ) mkl_jit_destroy( MatrixVectorProductJitterAlphaMinusOne ); + if ( MatrixVectorProductTranspJitterBetaOne != NULL ) mkl_jit_destroy( MatrixVectorProductTranspJitterBetaOne ); + if ( mkl_ipiv != NULL ) delete [] mkl_ipiv; #endif } @@ -231,22 +196,23 @@ void CSysMatrix::Initialize(unsigned long nPoint, unsigned long nPoi /*--- Generate MKL Kernels ---*/ -#if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) - /*--- Create MKL JIT kernels if not using adjoint solvers ---*/ - if (!config->GetContinuous_Adjoint() && !config->GetDiscrete_Adjoint()) - { - useMKL = true; +#ifdef USE_MKL + mkl_jit_create_dgemm( &MatrixMatrixProductJitter, MKL_ROW_MAJOR, MKL_NOTRANS, MKL_NOTRANS, nVar, nVar, nVar, 1.0, nVar, nVar, 0.0, nVar ); + MatrixMatrixProductKernel = mkl_jit_get_dgemm_ptr( MatrixMatrixProductJitter ); - mkl_jit_create_dgemm( &MatrixMatrixProductJitter, MKL_ROW_MAJOR, MKL_NOTRANS, MKL_NOTRANS, nVar, nVar, nVar, 1.0, nVar, nVar, 0.0, nVar ); - MatrixMatrixProductKernel = mkl_jit_get_dgemm_ptr( MatrixMatrixProductJitter ); + mkl_jit_create_dgemm( &MatrixVectorProductJitterBetaZero, MKL_COL_MAJOR, MKL_NOTRANS, MKL_NOTRANS, 1, nVar, nVar, 1.0, 1, nVar, 0.0, 1 ); + MatrixVectorProductKernelBetaZero = mkl_jit_get_dgemm_ptr( MatrixVectorProductJitterBetaZero ); - mkl_jit_create_dgemm( &MatrixVectorProductJitterBetaZero, MKL_COL_MAJOR, MKL_NOTRANS, MKL_NOTRANS, 1, nVar, nVar, 1.0, 1, nVar, 0.0, 1 ); - MatrixVectorProductKernelBetaZero = mkl_jit_get_dgemm_ptr( MatrixVectorProductJitterBetaZero ); + mkl_jit_create_dgemm( &MatrixVectorProductJitterBetaOne, MKL_COL_MAJOR, MKL_NOTRANS, MKL_NOTRANS, 1, nVar, nVar, 1.0, 1, nVar, 1.0, 1 ); + MatrixVectorProductKernelBetaOne = mkl_jit_get_dgemm_ptr( MatrixVectorProductJitterBetaOne ); - mkl_jit_create_dgemm( &MatrixVectorProductJitterBetaOne, MKL_COL_MAJOR, MKL_NOTRANS, MKL_NOTRANS, 1, nVar, nVar, 1.0, 1, nVar, 1.0, 1 ); - MatrixVectorProductKernelBetaOne = mkl_jit_get_dgemm_ptr( MatrixVectorProductJitterBetaOne ); - } + mkl_jit_create_dgemm( &MatrixVectorProductJitterAlphaMinusOne, MKL_COL_MAJOR, MKL_NOTRANS, MKL_NOTRANS, 1, nVar, nVar, -1.0, 1, nVar, 1.0, 1 ); + MatrixVectorProductKernelAlphaMinusOne = mkl_jit_get_dgemm_ptr( MatrixVectorProductJitterAlphaMinusOne ); + mkl_jit_create_dgemm( &MatrixVectorProductTranspJitterBetaOne, MKL_COL_MAJOR, MKL_NOTRANS, MKL_NOTRANS, nVar, 1, nVar, 1.0, nVar, nVar, 1.0, nVar ); + MatrixVectorProductTranspKernelBetaOne = mkl_jit_get_dgemm_ptr( MatrixVectorProductTranspJitterBetaOne ); + + mkl_ipiv = new lapack_int [ nVar ]; #endif /*--- Initialization matrix to zero ---*/ @@ -358,7 +324,6 @@ void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long block_weight = new ScalarType [nVar*nEqn]; block_inverse = new ScalarType [nVar*nEqn]; - prod_block_vector = new ScalarType [nEqn]; prod_row_vector = new ScalarType [nVar]; aux_vector = new ScalarType [nVar]; sum_vector = new ScalarType [nVar]; @@ -370,7 +335,6 @@ void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long for (iVar = 0; iVar < nVar*nEqn; iVar++) block_weight[iVar] = 0.0; for (iVar = 0; iVar < nVar*nEqn; iVar++) block_inverse[iVar] = 0.0; - for (iVar = 0; iVar < nEqn; iVar++) prod_block_vector[iVar] = 0.0; for (iVar = 0; iVar < nVar; iVar++) prod_row_vector[iVar] = 0.0; for (iVar = 0; iVar < nVar; iVar++) aux_vector[iVar] = 0.0; for (iVar = 0; iVar < nVar; iVar++) sum_vector[iVar] = 0.0; @@ -382,7 +346,6 @@ void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long if ((config->GetKind_Linear_Solver_Prec() == ILU) || ((config->GetKind_SU2() == SU2_DEF) && (config->GetKind_Deform_Linear_Solver_Prec() == ILU)) || ((config->GetKind_SU2() == SU2_DOT) && (config->GetKind_Deform_Linear_Solver_Prec() == ILU)) || - (config->GetKind_Linear_Solver() == SMOOTHER_ILU) || (config->GetFSI_Simulation() && config->GetKind_Deform_Linear_Solver_Prec() == ILU) || (config->GetDiscrete_Adjoint() && config->GetKind_DiscAdj_Linear_Prec() == ILU)) { @@ -391,6 +354,9 @@ void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long ILU_matrix = new ScalarType [nnz_ilu*nVar*nEqn]; for (iVar = 0; iVar < nnz_ilu*nVar*nEqn; iVar++) ILU_matrix[iVar] = 0.0; + invM = new ScalarType [nPointDomain*nVar*nEqn]; + for (iVar = 0; iVar < nPointDomain*nVar*nEqn; iVar++) invM[iVar] = 0.0; + } } @@ -401,15 +367,13 @@ void CSysMatrix::SetIndexes(unsigned long val_nPoint, unsigned long (config->GetKind_Linear_Solver_Prec() == LINELET) || ((config->GetKind_SU2() == SU2_DEF) && (config->GetKind_Deform_Linear_Solver_Prec() == JACOBI)) || ((config->GetKind_SU2() == SU2_DOT) && (config->GetKind_Deform_Linear_Solver_Prec() == JACOBI)) || - (config->GetKind_Linear_Solver() == SMOOTHER_JACOBI) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET) || (config->GetDiscrete_Adjoint() && config->GetKind_DiscAdj_Linear_Solver() == JACOBI) || (config->GetFSI_Simulation() && config->GetKind_Deform_Linear_Solver_Prec() == JACOBI)) { /*--- Reserve memory for the values of the inverse of the preconditioner. ---*/ - invM = new ScalarType [nPoint*nVar*nEqn]; - for (iVar = 0; iVar < nPoint*nVar*nEqn; iVar++) invM[iVar] = 0.0; + invM = new ScalarType [nPointDomain*nVar*nEqn]; + for (iVar = 0; iVar < nPointDomain*nVar*nEqn; iVar++) invM[iVar] = 0.0; } @@ -692,143 +656,6 @@ void CSysMatrix::CompleteComms(CSysVector & x, } -template -ScalarType *CSysMatrix::GetBlock(unsigned long block_i, unsigned long block_j) { - - unsigned long step = 0, index; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { return &(matrix[(row_ptr[block_i]+step-1)*nVar*nEqn]); } - } - return NULL; - -} - -template -ScalarType CSysMatrix::GetBlock(unsigned long block_i, unsigned long block_j, unsigned short iVar, unsigned short jVar) { - - unsigned long step = 0, index; - - for (index = row_ptr[block_i]; index < row_ptr[block_i+1]; index++) { - step++; - if (col_ind[index] == block_j) { return matrix[(row_ptr[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar]; } - } - return 0; - -} - -template -ScalarType *CSysMatrix::GetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j) { - - unsigned long step = 0, index; - - for (index = row_ptr_ilu[block_i]; index < row_ptr_ilu[block_i+1]; index++) { - step++; - if (col_ind_ilu[index] == block_j) { return &(ILU_matrix[(row_ptr_ilu[block_i]+step-1)*nVar*nEqn]); } - } - return NULL; - -} - -template -void CSysMatrix::SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block) { - - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr_ilu[block_i]; index < row_ptr_ilu[block_i+1]; index++) { - step++; - if (col_ind_ilu[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) - ILU_matrix[(row_ptr_ilu[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = val_block[iVar*nVar+jVar]; - break; - } - } - -} - -template -void CSysMatrix::SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block) { - - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr_ilu[block_i]; index < row_ptr_ilu[block_i+1]; index++) { - step++; - if (col_ind_ilu[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) - ILU_matrix[(row_ptr_ilu[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] = val_block[jVar*nVar+iVar]; - break; - } - } - -} - -template -void CSysMatrix::SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, ScalarType *val_block) { - - unsigned long iVar, jVar, index, step = 0; - - for (index = row_ptr_ilu[block_i]; index < row_ptr_ilu[block_i+1]; index++) { - step++; - if (col_ind_ilu[index] == block_j) { - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nEqn; jVar++) - ILU_matrix[(row_ptr_ilu[block_i]+step-1)*nVar*nEqn+iVar*nEqn+jVar] -= val_block[iVar*nVar+jVar]; - break; - } - } - -} - -template -void CSysMatrix::MatrixVectorProduct(ScalarType *matrix, ScalarType *vector, ScalarType *product) { - -#if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) - // NOTE: matrix/vector swapped due to column major kernel -- manual "CBLAS" setup. - if (useMKL) - { - MatrixVectorProductKernelBetaZero( MatrixVectorProductJitterBetaZero, vector, matrix, product ); - return; - } -#endif - - unsigned short iVar, jVar; - - for (iVar = 0; iVar < nVar; iVar++) { - product[iVar] = 0.0; - for (jVar = 0; jVar < nVar; jVar++) { - product[iVar] += matrix[iVar*nVar+jVar] * vector[jVar]; - } - } - -} - -template -void CSysMatrix::MatrixMatrixProduct(ScalarType *matrix_a, ScalarType *matrix_b, ScalarType *product) { - -#if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) - if (useMKL) - { - MatrixMatrixProductKernel( MatrixMatrixProductJitter, matrix_a, matrix_b, product ); - return; - } -#endif - - unsigned short iVar, jVar, kVar; - - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) { - product[iVar*nVar+jVar] = 0.0; - for (kVar = 0; kVar < nVar; kVar++) { - product[iVar*nVar+jVar] += matrix_a[iVar*nVar+kVar]*matrix_b[kVar*nVar+jVar]; - } - } - } - -} - template void CSysMatrix::DeleteValsRowi(unsigned long i) { @@ -846,324 +673,46 @@ void CSysMatrix::DeleteValsRowi(unsigned long i) { } template -ScalarType CSysMatrix::MatrixDeterminant(ScalarType **a, unsigned long n) { - - unsigned long i, j, j1, j2; - ScalarType det = 0; - ScalarType **m = NULL; - - if (n < 1) { } - else if (n == 1) { det = a[0][0]; } - else if (n == 2) { det = a[0][0] * a[1][1] - a[1][0] * a[0][1]; } - else { - det = 0.0; - - for (j1=0;j1 -void CSysMatrix::MatrixCoFactor(ScalarType **a, unsigned long n, ScalarType **b) { - - unsigned long i,j,ii,jj,i1,j1; - ScalarType det; - ScalarType **c; - - c = new ScalarType*[n-1]; - for (i=0;i -void CSysMatrix::MatrixTranspose(ScalarType **a, unsigned long n) { - - unsigned long i, j; - ScalarType tmp; - - for (i=1;i -void CSysMatrix::Gauss_Elimination(unsigned long block_i, ScalarType* rhs, bool transposed) { - - short iVar, jVar, kVar; // This is important, otherwise some compilers optimizations will fail - ScalarType weight, aux; - - ScalarType *Block = GetBlock(block_i, block_i); - - /*--- Copy block matrix, note that the original matrix - is modified by the algorithm---*/ - - if (!transposed) { - for (iVar = 0; iVar < (short)nVar; iVar++) - for (jVar = 0; jVar < (short)nVar; jVar++) - block[iVar*nVar+jVar] = Block[iVar*nVar+jVar]; - } else { - for (iVar = 0; iVar < (short)nVar; iVar++) - for (jVar = 0; jVar < (short)nVar; jVar++) - block[iVar*nVar+jVar] = Block[jVar*nVar+iVar]; - } - /*--- Gauss elimination ---*/ - - if (nVar == 1) { - rhs[0] /= block[0]; - } - else { - - /*--- Transform system in Upper Matrix ---*/ - - for (iVar = 1; iVar < (short)nVar; iVar++) { - for (jVar = 0; jVar < iVar; jVar++) { - weight = block[iVar*nVar+jVar] / block[jVar*nVar+jVar]; - for (kVar = jVar; kVar < (short)nVar; kVar++) - block[iVar*nVar+kVar] -= weight*block[jVar*nVar+kVar]; - rhs[iVar] -= weight*rhs[jVar]; - } - } - - /*--- Backwards substitution ---*/ - - rhs[nVar-1] = rhs[nVar-1] / block[nVar*nVar-1]; - for (iVar = (short)nVar-2; iVar >= 0; iVar--) { - aux = 0.0; - for (jVar = iVar+1; jVar < (short)nVar; jVar++) - aux += block[iVar*nVar+jVar]*rhs[jVar]; - rhs[iVar] = (rhs[iVar]-aux) / block[iVar*nVar+iVar]; - if (iVar == 0) break; - } - } - -} - -template -void CSysMatrix::Gauss_Elimination_ILUMatrix(unsigned long block_i, ScalarType* rhs) { - - short iVar, jVar, kVar; // This is important, otherwise some compilers optimizations will fail - ScalarType weight, aux; - - ScalarType *Block = GetBlock_ILUMatrix(block_i, block_i); - - /*--- Copy block matrix, note that the original matrix - is modified by the algorithm---*/ - - - // If source and dest overlap higher level problems occur, so memcpy is safe. And it is faster. - memcpy( block, Block, (nVar * nVar * sizeof(ScalarType)) ); - - //for (iVar = 0; iVar < (short)nVar; iVar++) - // for (jVar = 0; jVar < (short)nVar; jVar++) - // block[iVar*nVar+jVar] = Block[iVar*nVar+jVar]; - - /*--- Gauss elimination ---*/ - - if (nVar == 1) { - rhs[0] /= block[0]; - } - else { - -#if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) - if (useMKL) { - // With MKL_DIRECT_CALL enabled, this is significantly faster than native code on Intel Architectures. - lapack_int * ipiv = new lapack_int [ nVar ]; - LAPACKE_dgetrf( LAPACK_ROW_MAJOR, nVar, nVar, (double *)&block[0], nVar, ipiv ); - LAPACKE_dgetrs( LAPACK_ROW_MAJOR, 'N', nVar, 1, (double *)&block[0], nVar, ipiv, rhs, 1 ); - - delete [] ipiv; - return; - } -#endif - - /*--- Transform system in Upper Matrix ---*/ - - for (iVar = 1; iVar < (short)nVar; iVar++) { - for (jVar = 0; jVar < iVar; jVar++) { - weight = block[iVar*nVar+jVar] / block[jVar*nVar+jVar]; - for (kVar = jVar; kVar < (short)nVar; kVar++) - block[iVar*nVar+kVar] -= weight*block[jVar*nVar+kVar]; - rhs[iVar] -= weight*rhs[jVar]; - } - } - - /*--- Backwards substitution ---*/ - - rhs[nVar-1] = rhs[nVar-1] / block[nVar*nVar-1]; - for (iVar = (short)nVar-2; iVar >= 0; iVar--) { - aux = 0.0; - for (jVar = iVar+1; jVar < (short)nVar; jVar++) - aux += block[iVar*nVar+jVar]*rhs[jVar]; - rhs[iVar] = (rhs[iVar]-aux) / block[iVar*nVar+iVar]; - if (iVar == 0) break; - } - } - -} - -template -void CSysMatrix::Gauss_Elimination(ScalarType* Block, ScalarType* rhs) { - - short iVar, jVar, kVar; // This is important, otherwise some compilers optimizations will fail - ScalarType weight, aux; - - /*--- Copy block matrix, note that the original matrix - is modified by the algorithm---*/ - - for (iVar = 0; iVar < (short)nVar; iVar++) - for (jVar = 0; jVar < (short)nVar; jVar++) - block[iVar*nVar+jVar] = Block[iVar*nVar+jVar]; - - - if (nVar == 1) { - rhs[0] /= block[0]; - } - else { - /*--- Transform system in Upper Matrix ---*/ - for (iVar = 1; iVar < (short)nVar; iVar++) { - for (jVar = 0; jVar < iVar; jVar++) { - weight = block[iVar*nVar+jVar] / block[jVar*nVar+jVar]; - for (kVar = jVar; kVar < (short)nVar; kVar++) - block[iVar*nVar+kVar] -= weight*block[jVar*nVar+kVar]; - rhs[iVar] -= weight*rhs[jVar]; - } - } - - /*--- Backwards substitution ---*/ - rhs[nVar-1] = rhs[nVar-1] / block[nVar*nVar-1]; - for (iVar = (short)nVar-2; iVar >= 0; iVar--) { - aux = 0.0; - for (jVar = iVar+1; jVar < (short)nVar; jVar++) - aux += block[iVar*nVar+jVar]*rhs[jVar]; - rhs[iVar] = (rhs[iVar]-aux) / block[iVar*nVar+iVar]; - if (iVar == 0) break; - } - } - -} - -template -void CSysMatrix::ProdBlockVector(unsigned long block_i, unsigned long block_j, const CSysVector & vec) { - - unsigned long j = block_j*nVar; - unsigned short iVar, jVar; - - ScalarType *block = GetBlock(block_i, block_j); - - for (iVar = 0; iVar < nVar; iVar++) { - prod_block_vector[iVar] = 0; - for (jVar = 0; jVar < nVar; jVar++) - prod_block_vector[iVar] += block[iVar*nVar+jVar]*vec[j+jVar]; - } - -} - -template -void CSysMatrix::UpperProduct(CSysVector & vec, unsigned long row_i) { +void CSysMatrix::UpperProduct(const CSysVector & vec, unsigned long row_i) { - unsigned long iVar, index; + unsigned long iVar, index, col_j; for (iVar = 0; iVar < nVar; iVar++) prod_row_vector[iVar] = 0; for (index = row_ptr[row_i]; index < row_ptr[row_i+1]; index++) { - if (col_ind[index] > row_i) { - ProdBlockVector(row_i, col_ind[index], vec); - for (iVar = 0; iVar < nVar; iVar++) - prod_row_vector[iVar] += prod_block_vector[iVar]; + col_j = col_ind[index]; + if (col_j > row_i) { + MatrixVectorProductAdd(&matrix[index*nVar*nVar], &vec[col_j*nVar], prod_row_vector); } } } template -void CSysMatrix::LowerProduct(CSysVector & vec, unsigned long row_i) { +void CSysMatrix::LowerProduct(const CSysVector & vec, unsigned long row_i) { - unsigned long iVar, index; + unsigned long iVar, index, col_j; for (iVar = 0; iVar < nVar; iVar++) prod_row_vector[iVar] = 0; for (index = row_ptr[row_i]; index < row_ptr[row_i+1]; index++) { - if (col_ind[index] < row_i) { - ProdBlockVector(row_i, col_ind[index], vec); - for (iVar = 0; iVar < nVar; iVar++) - prod_row_vector[iVar] += prod_block_vector[iVar]; + col_j = col_ind[index]; + if (col_j < row_i) { + MatrixVectorProductAdd(&matrix[index*nVar*nVar], &vec[col_j*nVar], prod_row_vector); } } } template -void CSysMatrix::DiagonalProduct(CSysVector & vec, unsigned long row_i) { - - unsigned long iVar, index; - - for (iVar = 0; iVar < nVar; iVar++) - prod_row_vector[iVar] = 0; +void CSysMatrix::DiagonalProduct(const CSysVector & vec, unsigned long row_i) { - for (index = row_ptr[row_i]; index < row_ptr[row_i+1]; index++) { + for (unsigned long index = row_ptr[row_i]; index < row_ptr[row_i+1]; index++) { if (col_ind[index] == row_i) { - ProdBlockVector(row_i, col_ind[index], vec); - for (iVar = 0; iVar < nVar; iVar++) - prod_row_vector[iVar] += prod_block_vector[iVar]; + MatrixVectorProduct(&matrix[index*nVar*nVar], &vec[row_i*nVar], prod_row_vector); + break; } } @@ -1172,28 +721,14 @@ void CSysMatrix::DiagonalProduct(CSysVector & vec, unsig template void CSysMatrix::RowProduct(const CSysVector & vec, unsigned long row_i) { - unsigned long iVar, index; + unsigned long iVar, index, col_j; for (iVar = 0; iVar < nVar; iVar++) prod_row_vector[iVar] = 0; for (index = row_ptr[row_i]; index < row_ptr[row_i+1]; index++) { - ProdBlockVector(row_i, col_ind[index], vec); - for (iVar = 0; iVar < nVar; iVar++) - prod_row_vector[iVar] += prod_block_vector[iVar]; - } - -} - -template -void CSysMatrix::MatrixVectorProduct(const CSysVector & vec, CSysVector & prod) { - - unsigned long iPoint, iVar; - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - RowProduct(vec, iPoint); - for (iVar = 0; iVar < nVar; iVar++) - prod[iPoint*nVar+iVar] = prod_row_vector[iVar]; + col_j = col_ind[index]; + MatrixVectorProductAdd(&matrix[index*nVar*nVar], &vec[col_j*nVar], prod_row_vector); } } @@ -1201,7 +736,7 @@ void CSysMatrix::MatrixVectorProduct(const CSysVector & template void CSysMatrix::MatrixVectorProduct(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { - unsigned long prod_begin, vec_begin, mat_begin, index, iVar, jVar, row_i; + unsigned long prod_begin, vec_begin, mat_begin, index, row_i; /*--- Some checks for consistency between CSysMatrix and the CSysVectors ---*/ if ( (nVar != vec.GetNVar()) || (nVar != prod.GetNVar()) ) { @@ -1221,18 +756,7 @@ void CSysMatrix::MatrixVectorProduct(const CSysVector & for (index = row_ptr[row_i]; index < row_ptr[row_i+1]; index++) { vec_begin = col_ind[index]*nVar; // offset to beginning of block col_ind[index] mat_begin = (index*nVar*nVar); // offset to beginning of matrix block[row_i][col_ind[indx]] -#if defined(HAVE_MKL) && !(defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)) - if (useMKL) - { - MatrixVectorProductKernelBetaOne( MatrixVectorProductJitterBetaOne, (double *)&vec[ vec_begin ], (double *)&matrix[ mat_begin ], (double *)&prod[ prod_begin ] ); - continue; - } -#endif - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) { - prod[(unsigned long)(prod_begin+iVar)] += matrix[(unsigned long)(mat_begin+iVar*nVar+jVar)]*vec[(unsigned long)(vec_begin+jVar)]; - } - } + MatrixVectorProductAdd(&matrix[mat_begin], &vec[vec_begin], &prod[prod_begin]); } } @@ -1246,7 +770,7 @@ void CSysMatrix::MatrixVectorProduct(const CSysVector & template void CSysMatrix::MatrixVectorProductTransposed(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { - unsigned long prod_begin, vec_begin, mat_begin, index, iVar, jVar , row_i; + unsigned long prod_begin, vec_begin, mat_begin, index, row_i; /*--- Some checks for consistency between CSysMatrix and the CSysVectors ---*/ if ( (nVar != vec.GetNVar()) || (nVar != prod.GetNVar()) ) { @@ -1262,11 +786,7 @@ void CSysMatrix::MatrixVectorProductTransposed(const CSysVector::MatrixVectorProductTransposed(const CSysVector -void CSysMatrix::GetMultBlockBlock(ScalarType *c, ScalarType *a, ScalarType *b) { - - unsigned long iVar, jVar, kVar; - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) { - c[iVar*nVar+jVar] = 0.0; - for (kVar = 0; kVar < nVar; kVar++) - c[iVar*nVar+jVar] += a[iVar*nVar+kVar] * b[kVar*nVar+jVar]; - } - -} +void CSysMatrix::BuildJacobiPreconditioner(bool transpose) { -template -void CSysMatrix::GetMultBlockVector(ScalarType *c, ScalarType *a, ScalarType *b) { - - unsigned long iVar, jVar; - - for (iVar = 0; iVar < nVar; iVar++) { - c[iVar] = 0.0; - for (jVar = 0; jVar < nVar; jVar++) - c[iVar] += a[iVar*nVar+jVar] * b[jVar]; - } - -} + /*--- Build Jacobi preconditioner (M = D), compute and store the inverses of the diagonal blocks. ---*/ + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) + InverseDiagonalBlock(iPoint, &(invM[iPoint*nVar*nVar]), transpose); -template -void CSysMatrix::GetSubsBlock(ScalarType *c, ScalarType *a, ScalarType *b) { - - unsigned long iVar, jVar; - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - c[iVar*nVar+jVar] = a[iVar*nVar+jVar] - b[iVar*nVar+jVar]; - } template -void CSysMatrix::GetSubsVector(ScalarType *c, ScalarType *a, ScalarType *b) { - - unsigned long iVar; - - for (iVar = 0; iVar < nVar; iVar++) - c[iVar] = a[iVar] - b[iVar]; - -} +void CSysMatrix::ComputeJacobiPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { + + /*--- Apply Jacobi preconditioner, y = D^{-1} * x, the inverse of the diagonal is already known. ---*/ + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) + MatrixVectorProduct(&(invM[iPoint*nVar*nVar]), &vec[iPoint*nVar], &prod[iPoint*nVar]); + + /*--- MPI Parallelization ---*/ + InitiateComms(prod, geometry, config, SOLUTION_MATRIX); + CompleteComms(prod, geometry, config, SOLUTION_MATRIX); -template -void CSysMatrix::InverseBlock(ScalarType *Block, ScalarType *invBlock) { - - unsigned long iVar, jVar; - - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) - aux_vector[jVar] = 0.0; - aux_vector[iVar] = 1.0; - - /*--- Compute the i-th column of the inverse matrix ---*/ - Gauss_Elimination(Block, aux_vector); - - for (jVar = 0; jVar < nVar; jVar++) - invBlock[jVar*nVar+iVar] = aux_vector[jVar]; - } - } template -void CSysMatrix::InverseDiagonalBlock(unsigned long block_i, ScalarType *invBlock, bool transpose) { - - unsigned long iVar, jVar; - - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) - aux_vector[jVar] = 0.0; - aux_vector[iVar] = 1.0; - - /*--- Compute the i-th column of the inverse matrix ---*/ - - Gauss_Elimination(block_i, aux_vector, transpose); - for (jVar = 0; jVar < nVar; jVar++) - invBlock[jVar*nVar+iVar] = aux_vector[jVar]; - } - - // ScalarType Det, **Matrix, **CoFactor; - // ScalarType *Block = GetBlock(block_i, block_i); - // - // Matrix = new ScalarType*[nVar]; - // CoFactor = new ScalarType*[nVar]; - // for (iVar=0;iVar -void CSysMatrix::InverseDiagonalBlock_ILUMatrix(unsigned long block_i, ScalarType *invBlock) { - - unsigned long iVar, jVar; - - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) - aux_vector[jVar] = 0.0; - aux_vector[iVar] = 1.0; - - /*--- Compute the i-th column of the inverse matrix ---*/ - - Gauss_Elimination_ILUMatrix(block_i, aux_vector); - for (jVar = 0; jVar < nVar; jVar++) - invBlock[jVar*nVar+iVar] = aux_vector[jVar]; - } - - // ScalarType Det, **Matrix, **CoFactor; - // ScalarType *Block = GetBlock_ILUMatrix(block_i, block_i); - // - // Matrix = new ScalarType*[nVar]; - // CoFactor = new ScalarType*[nVar]; - // for (iVar=0;iVar -void CSysMatrix::BuildJacobiPreconditioner(bool transpose) { - - unsigned long iPoint, iVar, jVar; - - /*--- Compute Jacobi Preconditioner ---*/ - for (iPoint = 0; iPoint < nPoint; iPoint++) { - - /*--- Compute the inverse of the diagonal block ---*/ - InverseDiagonalBlock(iPoint, block_inverse, transpose); - - /*--- Set the inverse of the matrix to the invM structure (which is a vector) ---*/ - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - invM[iPoint*nVar*nVar+iVar*nVar+jVar] = block_inverse[iVar*nVar+jVar]; - } - -} - -template -void CSysMatrix::ComputeJacobiPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { - - unsigned long iPoint, iVar, jVar; - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - prod[(unsigned long)(iPoint*nVar+iVar)] = 0.0; - for (jVar = 0; jVar < nVar; jVar++) - prod[(unsigned long)(iPoint*nVar+iVar)] += - invM[(unsigned long)(iPoint*nVar*nVar+iVar*nVar+jVar)]*vec[(unsigned long)(iPoint*nVar+jVar)]; - } - } - - /*--- MPI Parallelization ---*/ - - InitiateComms(prod, geometry, config, SOLUTION_MATRIX); - CompleteComms(prod, geometry, config, SOLUTION_MATRIX); - -} - -template -unsigned long CSysMatrix::Jacobi_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config) { - - unsigned long iPoint, iVar, jVar; - ScalarType norm_r = 0.0, norm0 = 0.0; - int i = 0; - - /*--- Check the number of iterations requested ---*/ - - if (m < 1) { - char buf[100]; - SPRINTF(buf, "Illegal value for smoothing iterations, m = %lu", m ); - SU2_MPI::Error(string(buf), CURRENT_FUNCTION); - } - - /*--- Create vectors to hold the residual and the Matrix-Vector product - of the Jacobian matrix with the current solution (x^k). These must be - stored in order to perform multiple iterations of the smoother. ---*/ - - CSysVector r(b); - CSysVector A_x(b); - - /*--- Calculate the initial residual, compute norm, and check - if system is already solved. Recall, r holds b initially. ---*/ - - mat_vec(x, A_x); - r -= A_x; - - /*--- Only compute the residuals in full communication mode. ---*/ - - if (config->GetComm_Level() == COMM_FULL) { - - norm_r = r.norm(); - norm0 = b.norm(); - if ( (norm_r < tol*norm0) || (norm_r < eps) ) { - if (rank == MASTER_NODE) cout << "CSysMatrix::Jacobi_Smoother(): system solved by initial guess." << endl; - return 0; - } - - /*--- Set the norm to the initial residual value ---*/ - - norm0 = norm_r; - - /*--- Output header information including initial residual ---*/ - - if ((monitoring) && (rank == MASTER_NODE)) { - cout << "\n# " << "Jacobi Smoother" << " residual history" << endl; - cout << "# Residual tolerance target = " << tol << endl; - cout << "# Initial residual norm = " << norm_r << endl; - cout << " " << i << " " << norm_r/norm0 << endl; - } - - } - - /*--- Loop over all smoothing iterations ---*/ - - for (i = 0; i < (int)m; i++) { - - /*--- Apply the Jacobi smoother, i.e., multiply by the inverse of the - diagonal matrix of A, which was built in the preprocessing phase. Note - that we are directly updating the solution (x^k+1) during the loop. ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) - x[(unsigned long)(iPoint*nVar+iVar)] += - invM[(unsigned long)(iPoint*nVar*nVar+iVar*nVar+jVar)]*r[(unsigned long)(iPoint*nVar+jVar)]; - } - } - - /*--- MPI Parallelization ---*/ - - InitiateComms(x, geometry, config, SOLUTION_MATRIX); - CompleteComms(x, geometry, config, SOLUTION_MATRIX); - - /*--- Update the residual (r^k+1 = b - A*x^k+1) with the new solution ---*/ - - r = b; - mat_vec(x, A_x); - r -= A_x; - - /*--- Only compute the residuals in full communication mode. ---*/ - - if (config->GetComm_Level() == COMM_FULL) { - - /*--- Check if solution has converged, else output the relative - residual if necessary. ---*/ - - norm_r = r.norm(); - if (norm_r < tol*norm0) break; - if (((monitoring) && (rank == MASTER_NODE)) && ((i+1) % 5 == 0)) - cout << " " << i << " " << norm_r/norm0 << endl; - - } - - } - - if ((monitoring) && (rank == MASTER_NODE) && (config->GetComm_Level() == COMM_FULL)) { - cout << "# Jacobi smoother final (true) residual:" << endl; - cout << "# Iteration = " << i << ": |res|/|res0| = " << norm_r/norm0 << ".\n" << endl; - } - - return (unsigned long) i; - -} - -template -void CSysMatrix::BuildILUPreconditioner(bool transposed) { +void CSysMatrix::BuildILUPreconditioner(bool transposed) { unsigned long index, index_, iVar; - ScalarType *Block_ij, *Block_jk; + ScalarType *Block_ij; + const ScalarType *Block_jk; long iPoint, jPoint, kPoint; - /*--- Copy block matrix, note that the original matrix is modified by the algorithm, so that we have the factorization stored @@ -1623,6 +850,10 @@ void CSysMatrix::BuildILUPreconditioner(bool transposed) { for (iPoint = 1; iPoint < (long)nPointDomain; iPoint++) { + /*--- Invert and store the previous diagonal block to later compute the weight. ---*/ + + InverseDiagonalBlock_ILUMatrix(iPoint-1, &invM[(iPoint-1)*nVar*nVar]); + /*--- For each row (unknown), loop over all entries in A on this row row_ptr_ilu[iPoint+1] will have the index for the first entry on the next row. ---*/ @@ -1635,14 +866,13 @@ void CSysMatrix::BuildILUPreconditioner(bool transposed) { /*--- Check that this column is in the lower triangular portion ---*/ - if ((jPoint < iPoint) && (jPoint < (long)nPointDomain)) { + if (jPoint < iPoint) { - /*--- If we're in the lower triangle, get the pointer to this block, - invert it, and then right multiply against the original block ---*/ + /*--- If we're in the lower triangle, multiply the block by + the inverse of the corresponding diagonal block. ---*/ - Block_ij = GetBlock_ILUMatrix(iPoint, jPoint); - InverseDiagonalBlock_ILUMatrix(jPoint, block_inverse); - MatrixMatrixProduct(Block_ij, block_inverse, block_weight); + Block_ij = &ILU_matrix[index*nVar*nEqn]; + MatrixMatrixProduct(Block_ij, &invM[jPoint*nVar*nVar], block_weight); /*--- block_weight holds Aij*inv(Ajj). Jump to the row for jPoint ---*/ @@ -1656,9 +886,9 @@ void CSysMatrix::BuildILUPreconditioner(bool transposed) { upper triangular part, then multiply and modify the matrix. Here, Aik' = Aik - Aij*inv(Ajj)*Ajk. ---*/ - if ((kPoint >= jPoint) && (jPoint < (long)nPointDomain)) { + if (kPoint > jPoint) { - Block_jk = GetBlock_ILUMatrix(jPoint, kPoint); + Block_jk = &ILU_matrix[index_*nVar*nEqn]; MatrixMatrixProduct(block_weight, Block_jk, block); SubtractBlock_ILUMatrix(iPoint, kPoint, block); @@ -1668,30 +898,28 @@ void CSysMatrix::BuildILUPreconditioner(bool transposed) { /*--- Lastly, store block_weight in the lower triangular part, which will be reused during the forward solve in the precon/smoother. ---*/ - SetBlock_ILUMatrix(iPoint, jPoint, block_weight); + for (iVar = 0; iVar < nVar*nEqn; ++iVar) + Block_ij[iVar] = block_weight[iVar]; } } } + InverseDiagonalBlock_ILUMatrix(nPointDomain-1, &invM[(nPointDomain-1)*nVar*nVar]); + } template void CSysMatrix::ComputeILUPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { - unsigned long index; - ScalarType *Block_ij; + unsigned long index, iVar; + const ScalarType *Block_ij; long iPoint, jPoint; - unsigned short iVar; - /*--- Copy block matrix, note that the original matrix - is modified by the algorithm---*/ + /*--- Copy vector to then work on prod in place ---*/ - for (iPoint = 0; iPoint < (long)nPointDomain; iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - prod[iPoint*nVar+iVar] = vec[iPoint*nVar+iVar]; - } - } + for (iPoint = 0; iPoint < long(nPointDomain*nVar); iPoint++) + prod[iPoint] = vec[iPoint]; /*--- Forward solve the system using the lower matrix entries that were computed and stored during the ILU preprocessing. Note @@ -1700,39 +928,29 @@ void CSysMatrix::ComputeILUPreconditioner(const CSysVector= 0; iPoint--) { - for (iVar = 0; iVar < nVar; iVar++) sum_vector[iVar] = 0.0; + + for (iPoint = nPointDomain-1; iPoint >= 0; iPoint--) { + + for (iVar = 0; iVar < nVar; iVar++) + sum_vector[iVar] = prod[iPoint*nVar+iVar]; + for (index = row_ptr_ilu[iPoint]; index < row_ptr_ilu[iPoint+1]; index++) { jPoint = col_ind_ilu[index]; if ((jPoint >= iPoint+1) && (jPoint < (long)nPointDomain)) { - Block_ij = GetBlock_ILUMatrix(iPoint, jPoint); - MatrixVectorProduct(Block_ij, &prod[jPoint*nVar], aux_vector); - for (iVar = 0; iVar < nVar; iVar++) sum_vector[iVar] += aux_vector[iVar]; + Block_ij = &ILU_matrix[index*nVar*nEqn]; + MatrixVectorProductSub(Block_ij, &prod[jPoint*nVar], sum_vector); } } - for (iVar = 0; iVar < nVar; iVar++) prod[iPoint*nVar+iVar] = (prod[iPoint*nVar+iVar]-sum_vector[iVar]); - InverseDiagonalBlock_ILUMatrix(iPoint, block_inverse); - MatrixVectorProduct(block_inverse, &prod[iPoint*nVar], aux_vector); - for (iVar = 0; iVar < nVar; iVar++) prod[iPoint*nVar+iVar] = aux_vector[iVar]; - if (iPoint == 0) break; + + MatrixVectorProduct(&invM[iPoint*nVar*nVar], sum_vector, &prod[iPoint*nVar]); } /*--- MPI Parallelization ---*/ @@ -1742,166 +960,6 @@ void CSysMatrix::ComputeILUPreconditioner(const CSysVector -unsigned long CSysMatrix::ILU_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config) { - - unsigned long index; - ScalarType *Block_ij, omega = 1.0; - long iPoint, jPoint; - unsigned short iVar; - ScalarType norm_r = 0.0, norm0 = 0.0; - int i = 0; - - /*--- Check the number of iterations requested ---*/ - - if (m < 1) { - char buf[100]; - SPRINTF(buf, "Illegal value for smoothing iterations, m = %lu", m ); - SU2_MPI::Error(string(buf), CURRENT_FUNCTION); - } - - /*--- Create vectors to hold the residual and the Matrix-Vector product - of the Jacobian matrix with the current solution (x^k). These must be - stored in order to perform multiple iterations of the smoother. ---*/ - - CSysVector r(b); - CSysVector A_x(b); - - /*--- Calculate the initial residual, compute norm, and check - if system is already solved. Recall, r holds b initially. ---*/ - - mat_vec(x, A_x); - r -= A_x; - - /*--- Only compute the residuals in full communication mode. ---*/ - - if (config->GetComm_Level() == COMM_FULL) { - - norm_r = r.norm(); - norm0 = b.norm(); - if ( (norm_r < tol*norm0) || (norm_r < eps) ) { - if (rank == MASTER_NODE) cout << "CSysMatrix::ILU_Smoother(): system solved by initial guess." << endl; - return 0; - } - - /*--- Set the norm to the initial residual value ---*/ - - norm0 = norm_r; - - /*--- Output header information including initial residual ---*/ - - if ((monitoring) && (rank == MASTER_NODE)) { - cout << "\n# " << "ILU Smoother" << " residual history" << endl; - cout << "# Residual tolerance target = " << tol << endl; - cout << "# Initial residual norm = " << norm_r << endl; - cout << " " << i << " " << norm_r/norm0 << endl; - } - - } - - /*--- Loop over all smoothing iterations ---*/ - - for (i = 0; i < (int)m; i++) { - - /*--- Forward solve the system using the lower matrix entries that - were computed and stored during the ILU preprocessing. Note - that we are overwriting the residual vector as we go. ---*/ - - for (iPoint = 1; iPoint < (long)nPointDomain; iPoint++) { - - /*--- For each row (unknown), loop over all entries in A on this row - row_ptr_ilu[iPoint+1] will have the index for the first entry on the next - row. ---*/ - - for (index = row_ptr_ilu[iPoint]; index < row_ptr_ilu[iPoint+1]; index++) { - - /*--- jPoint here is the column for each entry on this row ---*/ - - jPoint = col_ind_ilu[index]; - - /*--- Check that this column is in the lower triangular portion ---*/ - - if ((jPoint < iPoint) && (jPoint < (long)nPointDomain)) { - - /*--- Lastly, get Aij*inv(Ajj) from the lower triangular part, which - was calculated in the preprocessing, and apply to r. ---*/ - - Block_ij = GetBlock_ILUMatrix(iPoint, jPoint); - MatrixVectorProduct(Block_ij, &r[jPoint*nVar], aux_vector); - for (iVar = 0; iVar < nVar; iVar++) - r[iPoint*nVar+iVar] -= aux_vector[iVar]; - - } - } - } - - /*--- Backwards substitution (starts at the last row) ---*/ - - InverseDiagonalBlock_ILUMatrix((nPointDomain-1), block_inverse); - MatrixVectorProduct(block_inverse, &r[(nPointDomain-1)*nVar], aux_vector); - - for (iVar = 0; iVar < nVar; iVar++) - r[(nPointDomain-1)*nVar + iVar] = aux_vector[iVar]; - - for (iPoint = nPointDomain-2; iPoint >= 0; iPoint--) { - for (iVar = 0; iVar < nVar; iVar++) sum_vector[iVar] = 0.0; - for (index = row_ptr_ilu[iPoint]; index < row_ptr_ilu[iPoint+1]; index++) { - jPoint = col_ind_ilu[index]; - if ((jPoint >= iPoint+1) && (jPoint < (long)nPointDomain)) { - Block_ij = GetBlock_ILUMatrix(iPoint, jPoint); - MatrixVectorProduct(Block_ij, &r[jPoint*nVar], aux_vector); - for (iVar = 0; iVar < nVar; iVar++) sum_vector[iVar] += aux_vector[iVar]; - } - } - for (iVar = 0; iVar < nVar; iVar++) r[iPoint*nVar+iVar] = (r[iPoint*nVar+iVar]-sum_vector[iVar]); - InverseDiagonalBlock_ILUMatrix(iPoint, block_inverse); - MatrixVectorProduct(block_inverse, &r[iPoint*nVar], aux_vector); - for (iVar = 0; iVar < nVar; iVar++) r[iPoint*nVar+iVar] = aux_vector[iVar]; - if (iPoint == 0) break; - } - - /*--- Update solution (x^k+1 = x^k + w*M^-1*r^k) using the residual vector, - which holds the update after applying the ILU smoother, i.e., M^-1*r^k. - Omega is a relaxation factor that we have currently set to 1.0. ---*/ - - x.Plus_AX(omega, r); - - /*--- MPI Parallelization ---*/ - - InitiateComms(x, geometry, config, SOLUTION_MATRIX); - CompleteComms(x, geometry, config, SOLUTION_MATRIX); - - /*--- Update the residual (r^k+1 = b - A*x^k+1) with the new solution ---*/ - - r = b; - mat_vec(x, A_x); - r -= A_x; - - /*--- Only compute the residuals in full communication mode. ---*/ - - if (config->GetComm_Level() == COMM_FULL) { - - /*--- Check if solution has converged, else output the relative - residual if necessary. ---*/ - - norm_r = r.norm(); - if (norm_r < tol*norm0) break; - if (((monitoring) && (rank == MASTER_NODE)) && ((i+1) % 5 == 0)) - cout << " " << i << " " << norm_r/norm0 << endl; - - } - - } - - if ((monitoring) && (rank == MASTER_NODE) && (config->GetComm_Level() == COMM_FULL)) { - cout << "# ILU smoother final (true) residual:" << endl; - cout << "# Iteration = " << i << ": |res|/|res0| = " << norm_r/norm0 << ".\n" << endl; - } - - return (unsigned int) i; - -} - template void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { unsigned long iPoint, iVar; @@ -1909,12 +967,10 @@ void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector::ComputeLU_SGSPreconditioner(const CSysVector= 0; iPoint--) { - DiagonalProduct(prod, iPoint); // Compute D.x* + DiagonalProduct(prod, iPoint); // Compute D.x* for (iVar = 0; iVar < nVar; iVar++) - aux_vector[iVar] = prod_row_vector[iVar]; // Compute aux_vector = D.x* - UpperProduct(prod, iPoint); // Compute U.x_(n+1) + aux_vector[iVar] = prod_row_vector[iVar]; // Compute aux_vector = D.x* + UpperProduct(prod, iPoint); // Compute U.x_(n+1) for (iVar = 0; iVar < nVar; iVar++) - aux_vector[iVar] -= prod_row_vector[iVar]; // Compute aux_vector = D.x*-U.x_(n+1) - Gauss_Elimination(iPoint, aux_vector); // Solve D.x* = aux_vector - for (iVar = 0; iVar < nVar; iVar++) - prod[iPoint*nVar + iVar] = aux_vector[iVar]; // Assesing x_(1) = solution + prod[iPoint*nVar+iVar] = aux_vector[iVar] - prod_row_vector[iVar]; // Compute aux_vector = D.x*-U.x_(n+1) + Gauss_Elimination(iPoint, &prod[iPoint*nVar]); // Solve D.x* = aux_vector } /*--- MPI Parallelization ---*/ @@ -1943,157 +997,20 @@ void CSysMatrix::ComputeLU_SGSPreconditioner(const CSysVector -unsigned long CSysMatrix::LU_SGS_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, ScalarType tol, unsigned long m, ScalarType *residual, bool monitoring, CGeometry *geometry, CConfig *config) { - - unsigned long iPoint, iVar; - ScalarType omega = 1.0; - ScalarType norm_r = 0.0, norm0 = 0.0; - int i = 0; - - /*--- Check the number of iterations requested ---*/ - - if (m < 1) { - char buf[100]; - SPRINTF(buf, "Illegal value for smoothing iterations, m = %lu", m ); - SU2_MPI::Error(string(buf), CURRENT_FUNCTION); - } - - /*--- Create vectors to hold the residual and the Matrix-Vector product - of the Jacobian matrix with the current solution (x^k). These must be - stored in order to perform multiple iterations of the smoother. ---*/ - - CSysVector r(b); - CSysVector A_x(b); - CSysVector xStar(x); - - /*--- Calculate the initial residual, compute norm, and check - if system is already solved. Recall, r holds b initially. ---*/ - - mat_vec(x, A_x); - r -= A_x; - - /*--- Only compute the residuals in full communication mode. ---*/ - - if (config->GetComm_Level() == COMM_FULL) { - - norm_r = r.norm(); - norm0 = b.norm(); - if ( (norm_r < tol*norm0) || (norm_r < eps) ) { - if (rank == MASTER_NODE) cout << "CSysMatrix::LU_SGS_Smoother(): system solved by initial guess." << endl; - return 0; - } - - /*--- Set the norm to the initial initial residual value ---*/ - - norm0 = norm_r; - - /*--- Output header information including initial residual ---*/ - - if ((monitoring) && (rank == MASTER_NODE)) { - cout << "\n# " << "LU_SGS Smoother" << " residual history" << endl; - cout << "# Residual tolerance target = " << tol << endl; - cout << "# Initial residual norm = " << norm_r << endl; - cout << " " << i << " " << norm_r/norm0 << endl; - } - - } - - /*--- Loop over all smoothing iterations ---*/ - - for (i = 0; i < (int)m; i++) { - - /*--- First part of the symmetric iteration: (D+L).x* = b ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - LowerProduct(xStar, iPoint); // Compute L.x* - for (iVar = 0; iVar < nVar; iVar++) - aux_vector[iVar] = r[iPoint*nVar+iVar] - prod_row_vector[iVar]; // Compute aux_vector = b - L.x* - Gauss_Elimination(iPoint, aux_vector); // Solve D.x* = aux_vector - for (iVar = 0; iVar < nVar; iVar++) - xStar[iPoint*nVar+iVar] = aux_vector[iVar]; // Assesing x* = solution, stored in r - } - - /*--- MPI Parallelization ---*/ - - InitiateComms(xStar, geometry, config, SOLUTION_MATRIX); - CompleteComms(xStar, geometry, config, SOLUTION_MATRIX); - - /*--- Second part of the symmetric iteration: (D+U).x_(1) = D.x* ---*/ - - for (iPoint = nPointDomain-1; (int)iPoint >= 0; iPoint--) { - DiagonalProduct(xStar, iPoint); // Compute D.x* - for (iVar = 0; iVar < nVar; iVar++) - aux_vector[iVar] = prod_row_vector[iVar]; // Compute aux_vector = D.x* - UpperProduct(xStar, iPoint); // Compute U.x_(n+1) - for (iVar = 0; iVar < nVar; iVar++) - aux_vector[iVar] -= prod_row_vector[iVar]; // Compute aux_vector = D.x*-U.x_(n+1) - Gauss_Elimination(iPoint, aux_vector); // Solve D.x* = aux_vector - for (iVar = 0; iVar < nVar; iVar++) - xStar[iPoint*nVar+iVar] = aux_vector[iVar]; // Assesing x_(1) = solution - } - - /*--- Update solution (x^k+1 = x^k + w*M^-1*r^k) using the xStar vector, - which holds the update after applying the LU_SGS smoother, i.e., M^-1*r^k. - Omega is a relaxation factor that we have currently set to 1.0. ---*/ - - x.Plus_AX(omega, xStar); - - /*--- MPI Parallelization ---*/ - - InitiateComms(x, geometry, config, SOLUTION_MATRIX); - CompleteComms(x, geometry, config, SOLUTION_MATRIX); - - /*--- Update the residual (r^k+1 = b - A*x^k+1) with the new solution ---*/ - - r = b; - mat_vec(x, A_x); - r -= A_x; - xStar = x; - - /*--- Only compute the residuals in full communication mode. ---*/ - - if (config->GetComm_Level() == COMM_FULL) { - - /*--- Check if solution has converged, else output the relative - residual if necessary. ---*/ - - norm_r = r.norm(); - if (norm_r < tol*norm0) break; - if (((monitoring) && (rank == MASTER_NODE)) && ((i+1) % 5 == 0)) - cout << " " << i << " " << norm_r/norm0 << endl; - - } - - } - - if ((monitoring) && (rank == MASTER_NODE) && (config->GetComm_Level() == COMM_FULL)) { - cout << "# LU_SGS smoother final (true) residual:" << endl; - cout << "# Iteration = " << i << ": |res|/|res0| = " << norm_r/norm0 << ".\n" << endl; - } - - return (unsigned int) i; - -} - template unsigned short CSysMatrix::BuildLineletPreconditioner(CGeometry *geometry, CConfig *config) { - bool *check_Point, add_point; + bool add_point; unsigned long iEdge, iPoint, jPoint, index_Point, iLinelet, iVertex, next_Point, counter, iElem; - unsigned short iMarker, iNode, ExtraLines = 100, MeanPoints; + unsigned short iMarker, iNode, MeanPoints; su2double alpha = 0.9, weight, max_weight, *normal, area, volume_iPoint, volume_jPoint; - unsigned long Local_nPoints, Local_nLineLets, Global_nPoints, Global_nLineLets; + unsigned long Local_nPoints, Local_nLineLets, Global_nPoints, Global_nLineLets, max_nElem; /*--- Memory allocation --*/ - check_Point = new bool [geometry->GetnPoint()]; - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) - check_Point[iPoint] = true; + vector check_Point(nPoint,true); - LineletBool = new bool[geometry->GetnPoint()]; - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint ++) - LineletBool[iPoint] = false; + LineletBool.resize(nPoint,false); nLinelet = 0; for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { @@ -2111,16 +1028,18 @@ unsigned short CSysMatrix::BuildLineletPreconditioner(CGeometry *geo /*--- Basic initial allocation ---*/ - LineletPoint = new vector[nLinelet + ExtraLines]; + LineletPoint.resize(nLinelet); /*--- Define the basic linelets, starting from each vertex ---*/ + iLinelet = 0; + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { if ((config->GetMarker_All_KindBC(iMarker) == HEAT_FLUX ) || (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL ) || (config->GetMarker_All_KindBC(iMarker) == EULER_WALL ) || - (config->GetMarker_All_KindBC(iMarker) == DISPLACEMENT_BOUNDARY)) { - iLinelet = 0; + (config->GetMarker_All_KindBC(iMarker) == DISPLACEMENT_BOUNDARY)) + { for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); LineletPoint[iLinelet].push_back(iPoint); @@ -2242,28 +1161,9 @@ unsigned short CSysMatrix::BuildLineletPreconditioner(CGeometry *geo /*--- Memory allocation --*/ - UBlock = new ScalarType* [max_nElem]; - invUBlock = new ScalarType* [max_nElem]; - LBlock = new ScalarType* [max_nElem]; - yVector = new ScalarType* [max_nElem]; - zVector = new ScalarType* [max_nElem]; - rVector = new ScalarType* [max_nElem]; - for (iElem = 0; iElem < max_nElem; iElem++) { - UBlock[iElem] = new ScalarType [nVar*nVar]; - invUBlock[iElem] = new ScalarType [nVar*nVar]; - LBlock[iElem] = new ScalarType [nVar*nVar]; - yVector[iElem] = new ScalarType [nVar]; - zVector[iElem] = new ScalarType [nVar]; - rVector[iElem] = new ScalarType [nVar]; - } - - LFBlock = new ScalarType [nVar*nVar]; - LyVector = new ScalarType [nVar]; - FzVector = new ScalarType [nVar]; - - /*--- Memory deallocation --*/ - - delete [] check_Point; + LineletUpper.resize(max_nElem,NULL); + LineletInvDiag.resize(max_nElem*nVar*nVar,0.0); + LineletVector.resize(max_nElem*nVar,0.0); return MeanPoints; @@ -2272,133 +1172,163 @@ unsigned short CSysMatrix::BuildLineletPreconditioner(CGeometry *geo template void CSysMatrix::ComputeLineletPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config) { - - unsigned long iVar, jVar, nElem = 0, iLinelet, im1Point, iPoint, ip1Point, iElem; - long iElemLoop; - ScalarType *block; - - if (size == SINGLE_NODE) { - - /*--- Jacobi preconditioning if there is no linelet ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - if (!LineletBool[iPoint]) { - for (iVar = 0; iVar < nVar; iVar++) { - prod[(unsigned long)(iPoint*nVar+iVar)] = 0.0; - for (jVar = 0; jVar < nVar; jVar++) - prod[(unsigned long)(iPoint*nVar+iVar)] += - invM[(unsigned long)(iPoint*nVar*nVar+iVar*nVar+jVar)]*vec[(unsigned long)(iPoint*nVar+jVar)]; - } - } + + unsigned long iVar, iElem, nElem, iLinelet, iPoint, im1Point; + /*--- Pointers to lower, upper, and diagonal blocks ---*/ + const ScalarType *l = NULL, *u = NULL, *d = NULL; + /*--- Inverse of d_{i-1}, modified d_i, modified b_i (rhs) ---*/ + ScalarType *inv_dm1 = NULL, *d_prime = NULL, *b_prime = NULL; + +// if (size != SINGLE_NODE) +// SU2_MPI::Error("Linelet not implemented in parallel.", CURRENT_FUNCTION); + + /*--- Jacobi preconditioning where there is no linelet ---*/ + + for (iPoint = 0; iPoint < nPointDomain; iPoint++) + if (!LineletBool[iPoint]) + MatrixVectorProduct(&(invM[iPoint*nVar*nVar]), &vec[iPoint*nVar], &prod[iPoint*nVar]); + + /*--- MPI Parallelization ---*/ + + InitiateComms(prod, geometry, config, SOLUTION_MATRIX); + CompleteComms(prod, geometry, config, SOLUTION_MATRIX); + + /*--- Solve linelet using the Thomas algorithm ---*/ + + for (iLinelet = 0; iLinelet < nLinelet; iLinelet++) { + + nElem = LineletPoint[iLinelet].size(); + + /*--- Initialize the solution vector with the rhs ---*/ + + for (iElem = 0; iElem < nElem; iElem++) { + iPoint = LineletPoint[iLinelet][iElem]; + for (iVar = 0; iVar < nVar; iVar++) + LineletVector[iElem*nVar+iVar] = vec[iPoint*nVar+iVar]; } - - /*--- MPI Parallelization ---*/ - - InitiateComms(prod, geometry, config, SOLUTION_MATRIX); - CompleteComms(prod, geometry, config, SOLUTION_MATRIX); - - /*--- Solve linelet using a Thomas' algorithm ---*/ - - for (iLinelet = 0; iLinelet < nLinelet; iLinelet++) { - - nElem = LineletPoint[iLinelet].size(); - - /*--- Copy vec vector to the new structure ---*/ - - for (iElem = 0; iElem < nElem; iElem++) { - iPoint = LineletPoint[iLinelet][iElem]; - for (iVar = 0; iVar < nVar; iVar++) - rVector[iElem][iVar] = vec[(unsigned long)(iPoint*nVar+iVar)]; - } - - /*--- Initialization (iElem = 0) ---*/ - - iPoint = LineletPoint[iLinelet][0]; - block = GetBlock(iPoint, iPoint); - for (iVar = 0; iVar < nVar; iVar++) { - yVector[0][iVar] = rVector[0][iVar]; - for (jVar = 0; jVar < nVar; jVar++) - UBlock[0][iVar*nVar+jVar] = block[iVar*nVar+jVar]; - } - - /*--- Main loop (without iElem = 0) ---*/ - - for (iElem = 1; iElem < nElem; iElem++) { - - im1Point = LineletPoint[iLinelet][iElem-1]; - iPoint = LineletPoint[iLinelet][iElem]; - - InverseBlock(UBlock[iElem-1], invUBlock[iElem-1]); - block = GetBlock(iPoint, im1Point); GetMultBlockBlock(LBlock[iElem], block, invUBlock[iElem-1]); - block = GetBlock(im1Point, iPoint); GetMultBlockBlock(LFBlock, LBlock[iElem], block); - block = GetBlock(iPoint, iPoint); GetSubsBlock(UBlock[iElem], block, LFBlock); - - /*--- Forward substituton ---*/ - - GetMultBlockVector(LyVector, LBlock[iElem], yVector[iElem-1]); - GetSubsVector(yVector[iElem], rVector[iElem], LyVector); - - } - - /*--- Backward substituton ---*/ - - InverseBlock(UBlock[nElem-1], invUBlock[nElem-1]); - GetMultBlockVector(zVector[nElem-1], invUBlock[nElem-1], yVector[nElem-1]); - - for (iElemLoop = nElem-2; iElemLoop >= 0; iElemLoop--) { - iPoint = LineletPoint[iLinelet][iElemLoop]; - ip1Point = LineletPoint[iLinelet][iElemLoop+1]; - block = GetBlock(iPoint, ip1Point); GetMultBlockVector(FzVector, block, zVector[iElemLoop+1]); - GetSubsVector(aux_vector, yVector[iElemLoop], FzVector); - GetMultBlockVector(zVector[iElemLoop], invUBlock[iElemLoop], aux_vector); - } - - /*--- Copy zVector to the prod vector ---*/ - - for (iElem = 0; iElem < nElem; iElem++) { - iPoint = LineletPoint[iLinelet][iElem]; - for (iVar = 0; iVar < nVar; iVar++) - prod[(unsigned long)(iPoint*nVar+iVar)] = zVector[iElem][iVar]; - } - + + /*--- Forward pass, eliminate lower entries, modify diagonal and rhs ---*/ + + iPoint = LineletPoint[iLinelet][0]; + d = GetBlock(iPoint, iPoint); + for (iVar = 0; iVar < nVar*nVar; ++iVar) + LineletInvDiag[iVar] = d[iVar]; + + for (iElem = 1; iElem < nElem; iElem++) { + + /*--- Setup pointers to required matrices and vectors ---*/ + im1Point = LineletPoint[iLinelet][iElem-1]; + iPoint = LineletPoint[iLinelet][iElem]; + + d = GetBlock(iPoint, iPoint); + l = GetBlock(iPoint, im1Point); + u = GetBlock(im1Point, iPoint); + + inv_dm1 = &LineletInvDiag[(iElem-1)*nVar*nVar]; + d_prime = &LineletInvDiag[iElem*nVar*nVar]; + b_prime = &LineletVector[iElem*nVar]; + + /*--- Invert previous modified diagonal ---*/ + MatrixInverse(inv_dm1, inv_dm1); + + /*--- Left-multiply by lower block to obtain the weight ---*/ + MatrixMatrixProduct(l, inv_dm1, block_weight); + + /*--- Multiply weight by upper block to modify current diagonal ---*/ + MatrixMatrixProduct(block_weight, u, d_prime); + MatrixSubtraction(d, d_prime, d_prime); + + /*--- Update the rhs ---*/ + MatrixVectorProduct(block_weight, &LineletVector[(iElem-1)*nVar], aux_vector); + VectorSubtraction(b_prime, aux_vector, b_prime); + + /*--- Cache upper block pointer for the backward substitution phase ---*/ + LineletUpper[iElem-1] = u; } - - /*--- MPI Parallelization ---*/ - - InitiateComms(prod, geometry, config, SOLUTION_MATRIX); - CompleteComms(prod, geometry, config, SOLUTION_MATRIX); - - } - else { - SU2_MPI::Error("Linelet not implemented in parallel.", CURRENT_FUNCTION); + + /*--- Backwards substitution, LineletVector becomes the solution ---*/ + + /*--- x_n = d_n^{-1} * b_n ---*/ + Gauss_Elimination(&LineletInvDiag[(nElem-1)*nVar*nVar], &LineletVector[(nElem-1)*nVar]); + + /*--- x_i = d_i^{-1}*(b_i - u_i*x_{i+1}) ---*/ + for (iElem = nElem-1; iElem > 0; --iElem) { + inv_dm1 = &LineletInvDiag[(iElem-1)*nVar*nVar]; + MatrixVectorProduct(LineletUpper[iElem-1], &LineletVector[iElem*nVar], aux_vector); + VectorSubtraction(&LineletVector[(iElem-1)*nVar], aux_vector, aux_vector); + MatrixVectorProduct(inv_dm1, aux_vector, &LineletVector[(iElem-1)*nVar]); + } + + /*--- Copy results to product vector ---*/ + + for (iElem = 0; iElem < nElem; iElem++) { + iPoint = LineletPoint[iLinelet][iElem]; + for (iVar = 0; iVar < nVar; iVar++) + prod[iPoint*nVar+iVar] = LineletVector[iElem*nVar+iVar]; + } + } - + + /*--- MPI Parallelization ---*/ + + InitiateComms(prod, geometry, config, SOLUTION_MATRIX); + CompleteComms(prod, geometry, config, SOLUTION_MATRIX); + } template void CSysMatrix::ComputeResidual(const CSysVector & sol, const CSysVector & f, CSysVector & res) { - unsigned long iPoint, iVar; - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { RowProduct(sol, iPoint); - for (iVar = 0; iVar < nVar; iVar++) { - res[iPoint*nVar+iVar] = prod_row_vector[iVar] - f[iPoint*nVar+iVar]; - } + VectorSubtraction(prod_row_vector, &f[iPoint*nVar], &res[iPoint*nVar]); } } +template +template +void CSysMatrix::EnforceSolutionAtNode(const unsigned long node_i, const OtherType *x_i, CSysVector & b) { + + /*--- Both row and column associated with node i are eliminated (Block_ii = I and all else 0) to preserve eventual symmetry. ---*/ + /*--- The vector is updated with the product of column i by the known (enforced) solution at node i. ---*/ + + unsigned long iPoint, iVar, jVar, index, mat_begin; + + /*--- Delete whole row first. ---*/ + for (index = row_ptr[node_i]*nVar*nVar; index < row_ptr[node_i+1]*nVar*nVar; ++index) + matrix[index] = 0.0; + + /*--- Update b with the column product and delete column. ---*/ + for (iPoint = 0; iPoint < nPoint; ++iPoint) { + for (index = row_ptr[iPoint]; index < row_ptr[iPoint+1]; ++index) { + if (col_ind[index] == node_i) + { + mat_begin = index*nVar*nVar; + + for(iVar = 0; iVar < nVar; ++iVar) + for(jVar = 0; jVar < nVar; ++jVar) + b[iPoint*nVar+iVar] -= matrix[mat_begin+iVar*nVar+jVar] * x_i[jVar]; + + /*--- If on diagonal, set diagonal of block to 1, else delete block. ---*/ + if (iPoint == node_i) + for (iVar = 0; iVar < nVar; ++iVar) matrix[mat_begin+iVar*(nVar+1)] = 1.0; + else + for (iVar = 0; iVar < nVar*nVar; iVar++) matrix[mat_begin+iVar] = 0.0; + } + } + } + + /*--- Set know solution in rhs vector. ---*/ + for (iVar = 0; iVar < nVar; iVar++) b[node_i*nVar+iVar] = x_i[iVar]; + +} + /*--- Explicit instantiations ---*/ template class CSysMatrix; template void CSysMatrix::InitiateComms(CSysVector&, CGeometry*, CConfig*, unsigned short); template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short); -template class CSysMatrixVectorProduct; -template class CSysMatrixVectorProductTransposed; -template class CJacobiPreconditioner; -template class CILUPreconditioner; -template class CLU_SGSPreconditioner; -template class CLineletPreconditioner; +template void CSysMatrix::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector&); #ifdef CODI_REVERSE_TYPE template class CSysMatrix; @@ -2406,10 +1336,6 @@ template void CSysMatrix::InitiateComms(CSysVector::InitiateComms(CSysVector&, CGeometry*, CConfig*, unsigned short); template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short); template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short); -template class CSysMatrixVectorProduct; -template class CSysMatrixVectorProductTransposed; -template class CJacobiPreconditioner; -template class CILUPreconditioner; -template class CLU_SGSPreconditioner; -template class CLineletPreconditioner; +template void CSysMatrix::EnforceSolutionAtNode(unsigned long, const passivedouble*, CSysVector&); +template void CSysMatrix::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector&); #endif diff --git a/Common/src/vector_structure.cpp b/Common/src/vector_structure.cpp index 77ac968741f6..61a56418d199 100644 --- a/Common/src/vector_structure.cpp +++ b/Common/src/vector_structure.cpp @@ -423,6 +423,40 @@ ScalarType *CSysVector::GetBlock(unsigned long val_ipoint) { return &vec_val[val_ipoint*nVar]; } +template +template +void CSysVector::PassiveCopy(const CSysVector& other) { + + /*--- This is a method and not the overload of an operator to make sure who + calls it knows the consequence to the derivative information (lost) ---*/ + + /*--- check if self-assignment, otherwise perform deep copy ---*/ + if ((const void*)this == (const void*)&other) return; + + /*--- determine if (re-)allocation is needed ---*/ + if (nElm != other.GetLocSize() && vec_val != NULL) { + delete [] vec_val; + vec_val = NULL; + } + + /*--- copy ---*/ + nElm = other.GetLocSize(); + nElmDomain = other.GetNElmDomain(); + nBlk = other.GetNBlk(); + nBlkDomain = other.GetNBlkDomain(); + nVar = other.GetNVar(); + + if (vec_val == NULL) + vec_val = new ScalarType[nElm]; + + for (unsigned long i = 0; i < nElm; i++) + vec_val[i] = SU2_TYPE::GetValue(other[i]); + +#ifdef HAVE_MPI + nElmGlobal = other.GetSize(); +#endif +} + template ScalarType dotProd(const CSysVector & u, const CSysVector & v) { @@ -450,10 +484,13 @@ ScalarType dotProd(const CSysVector & u, const CSysVector; template CSysVector operator*(const su2double&, const CSysVector&); +template void CSysVector::PassiveCopy(const CSysVector&); template su2double dotProd(const CSysVector & u, const CSysVector & v); #ifdef CODI_REVERSE_TYPE template class CSysVector; template CSysVector operator*(const passivedouble&, const CSysVector&); +template void CSysVector::PassiveCopy(const CSysVector&); +template void CSysVector::PassiveCopy(const CSysVector&); template passivedouble dotProd(const CSysVector & u, const CSysVector & v); #endif diff --git a/SU2_CFD/src/solver_adjoint_mean.cpp b/SU2_CFD/src/solver_adjoint_mean.cpp index d7b63dfcedc8..78d556926c56 100644 --- a/SU2_CFD/src/solver_adjoint_mean.cpp +++ b/SU2_CFD/src/solver_adjoint_mean.cpp @@ -187,8 +187,7 @@ CAdjEulerSolver::CAdjEulerSolver(CGeometry *geometry, CConfig *config, unsigned cout << "Initialize Jacobian structure (Adjoint Euler). MG level: " << iMesh <<"." << endl; Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } @@ -5292,8 +5291,7 @@ CAdjNSSolver::CAdjNSSolver(CGeometry *geometry, CConfig *config, unsigned short cout << "Initialize Jacobian structure (Adjoint N-S). MG level: " << iMesh <<"." << endl; Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } diff --git a/SU2_CFD/src/solver_adjoint_turbulent.cpp b/SU2_CFD/src/solver_adjoint_turbulent.cpp index 1b58346e3b26..82fc48c85870 100644 --- a/SU2_CFD/src/solver_adjoint_turbulent.cpp +++ b/SU2_CFD/src/solver_adjoint_turbulent.cpp @@ -96,8 +96,7 @@ CAdjTurbSolver::CAdjTurbSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- Initialization of the structure of the whole Jacobian ---*/ Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } diff --git a/SU2_CFD/src/solver_direct_elasticity.cpp b/SU2_CFD/src/solver_direct_elasticity.cpp index f01c228332b0..2f2b033f8f58 100644 --- a/SU2_CFD/src/solver_direct_elasticity.cpp +++ b/SU2_CFD/src/solver_direct_elasticity.cpp @@ -1973,82 +1973,35 @@ void CFEASolver::BC_Clamped(CGeometry *geometry, CSolver **solver_container, CNu unsigned short val_marker) { unsigned long iPoint, iVertex; - unsigned long iVar, jVar; bool dynamic = (config->GetDynamic_Analysis() == DYNAMIC); + + Solution[0] = 0.0; + Solution[1] = 0.0; + if (nDim==3) Solution[2] = 0.0; for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Get node index ---*/ - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - if (geometry->node[iPoint]->GetDomain()) { - - if (nDim == 2) { - Solution[0] = 0.0; Solution[1] = 0.0; - Residual[0] = 0.0; Residual[1] = 0.0; - } - else { - Solution[0] = 0.0; Solution[1] = 0.0; Solution[2] = 0.0; - Residual[0] = 0.0; Residual[1] = 0.0; Residual[2] = 0.0; - } - - node[iPoint]->SetSolution(Solution); - - if (dynamic) { - node[iPoint]->SetSolution_Vel(Solution); - node[iPoint]->SetSolution_Accel(Solution); - } - - - /*--- Initialize the reaction vector ---*/ - LinSysReact.SetBlock(iPoint, Residual); - - LinSysRes.SetBlock(iPoint, Residual); - LinSysSol.SetBlock(iPoint, Solution); - - /*--- STRONG ENFORCEMENT OF THE DISPLACEMENT BOUNDARY CONDITION ---*/ - - /*--- Delete the columns for a particular node ---*/ - - for (iVar = 0; iVar < nPoint; iVar++) { - if (iVar==iPoint) { - Jacobian.SetBlock(iVar,iPoint,mId_Aux); - } - else { - Jacobian.SetBlock(iVar,iPoint,mZeros_Aux); - } - } - - /*--- Delete the rows for a particular node ---*/ - for (jVar = 0; jVar < nPoint; jVar++) { - if (iPoint!=jVar) { - Jacobian.SetBlock(iPoint,jVar,mZeros_Aux); - } - } - - /*--- If the problem is dynamic ---*/ - /*--- Enforce that in the previous time step all nodes had 0 U, U', U'' ---*/ - - if(dynamic) { - - node[iPoint]->SetSolution_time_n(Solution); - node[iPoint]->SetSolution_Vel_time_n(Solution); - node[iPoint]->SetSolution_Accel_time_n(Solution); - - } - - } else { - - /*--- Delete the column (iPoint is halo so Send/Recv does the rest) ---*/ - - for (iVar = 0; iVar < nPoint; iVar++) Jacobian.SetBlock(iVar,iPoint,mZeros_Aux); - + + /*--- Set and enforce solution at current and previous time-step ---*/ + node[iPoint]->SetSolution(Solution); + + if (dynamic) { + node[iPoint]->SetSolution_Vel(Solution); + node[iPoint]->SetSolution_Accel(Solution); + node[iPoint]->SetSolution_time_n(Solution); + node[iPoint]->SetSolution_Vel_time_n(Solution); + node[iPoint]->SetSolution_Accel_time_n(Solution); } - + + LinSysSol.SetBlock(iPoint, Solution); + LinSysReact.SetBlock(iPoint, Solution); + Jacobian.EnforceSolutionAtNode(iPoint, Solution, LinSysRes); + } - + } void CFEASolver::BC_Clamped_Post(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, CConfig *config, @@ -2084,8 +2037,8 @@ void CFEASolver::BC_Clamped_Post(CGeometry *geometry, CSolver **solver_container void CFEASolver::BC_DispDir(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, CConfig *config, unsigned short val_marker) { - unsigned short iDim, jDim; - + unsigned short iDim; + unsigned long iNode, iVertex; su2double DispDirVal = config->GetDisp_Dir_Value(config->GetMarker_All_TagBound(val_marker)); su2double DispDirMult = config->GetDisp_Dir_Multiplier(config->GetMarker_All_TagBound(val_marker)); su2double *Disp_Dir_Local= config->GetDisp_Dir(config->GetMarker_All_TagBound(val_marker)); @@ -2147,84 +2100,14 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, CSolver **solver_container, CNu for (iDim = 0; iDim < nDim; iDim++) Disp_Dir[iDim] = TotalDisp * Disp_Dir_Unit[iDim]; - unsigned long iNode, iVertex; - unsigned long iPoint, jPoint; - - su2double valJacobian_ij_00 = 0.0; - su2double auxJacobian_ij[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Get node index ---*/ - iNode = geometry->vertex[val_marker][iVertex]->GetNode(); - if (geometry->node[iNode]->GetDomain()) { - - if (nDim == 2) { - Solution[0] = Disp_Dir[0]; Solution[1] = Disp_Dir[1]; - Residual[0] = Disp_Dir[0]; Residual[1] = Disp_Dir[1]; - } - else { - Solution[0] = Disp_Dir[0]; Solution[1] = Disp_Dir[1]; Solution[2] = Disp_Dir[2]; - Residual[0] = Disp_Dir[0]; Residual[1] = Disp_Dir[1]; Residual[2] = Disp_Dir[2]; - } - - /*--- Initialize the reaction vector ---*/ - - LinSysRes.SetBlock(iNode, Residual); - LinSysSol.SetBlock(iNode, Solution); - - /*--- STRONG ENFORCEMENT OF THE DISPLACEMENT BOUNDARY CONDITION ---*/ - - /*--- Delete the full row for node iNode ---*/ - for (jPoint = 0; jPoint < nPoint; jPoint++){ - if (iNode != jPoint) { - Jacobian.SetBlock(iNode,jPoint,mZeros_Aux); - } - else{ - Jacobian.SetBlock(iNode,jPoint,mId_Aux); - } - } - - } - - /*--- Always delete the iNode column, even for halos ---*/ - - for (iPoint = 0; iPoint < nPoint; iPoint++) { - - /*--- Check if the term K(iPoint, iNode) is 0 ---*/ - valJacobian_ij_00 = Jacobian.GetBlock(iPoint,iNode,0,0); - - /*--- If the node iNode has a crossed dependency with the point iPoint ---*/ - if (valJacobian_ij_00 != 0.0 ){ - - /*--- Retrieve the Jacobian term ---*/ - for (iDim = 0; iDim < nDim; iDim++){ - for (jDim = 0; jDim < nDim; jDim++){ - auxJacobian_ij[iDim][jDim] = Jacobian.GetBlock(iPoint,iNode,iDim,jDim); - } - } - - /*--- Multiply by the imposed displacement ---*/ - for (iDim = 0; iDim < nDim; iDim++){ - Residual[iDim] = 0.0; - for (jDim = 0; jDim < nDim; jDim++){ - Residual[iDim] += auxJacobian_ij[iDim][jDim] * Disp_Dir[jDim]; - } - } - - if (iNode != iPoint) { - /*--- The term is substracted from the residual (right hand side) ---*/ - LinSysRes.SubtractBlock(iPoint, Residual); - /*--- The Jacobian term is now set to 0 ---*/ - Jacobian.SetBlock(iPoint,iNode,mZeros_Aux); - } - - } - - } - + /*--- Set and enforce solution ---*/ + LinSysSol.SetBlock(iNode, Disp_Dir); + Jacobian.EnforceSolutionAtNode(iNode, Disp_Dir, LinSysRes); } } diff --git a/SU2_CFD/src/solver_direct_heat.cpp b/SU2_CFD/src/solver_direct_heat.cpp index 917b468d7fd7..ed32ad7af6fc 100644 --- a/SU2_CFD/src/solver_direct_heat.cpp +++ b/SU2_CFD/src/solver_direct_heat.cpp @@ -130,8 +130,7 @@ CHeatSolverFVM::CHeatSolverFVM(CGeometry *geometry, CConfig *config, unsigned sh if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (heat equation) MG level: " << iMesh << "." << endl; Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index cc52d8ccbca2..9ceb6b6f83d0 100644 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -476,8 +476,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, unsigned short if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (Euler). MG level: " << iMesh <<"." << endl; Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } @@ -14776,8 +14775,7 @@ CNSSolver::CNSSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (Navier-Stokes). MG level: " << iMesh <<"." << endl; Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index b213f997db0e..9670cd473f58 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -324,8 +324,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (Euler). MG level: " << iMesh <<"." << endl; Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } @@ -6989,8 +6988,7 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (Navier-Stokes). MG level: " << iMesh <<"." << endl; Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } diff --git a/SU2_CFD/src/solver_direct_transition.cpp b/SU2_CFD/src/solver_direct_transition.cpp index 231dacf4eb6f..64349e88bf41 100644 --- a/SU2_CFD/src/solver_direct_transition.cpp +++ b/SU2_CFD/src/solver_direct_transition.cpp @@ -108,8 +108,7 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- Initialization of the structure of the whole Jacobian ---*/ Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } diff --git a/SU2_CFD/src/solver_direct_turbulent.cpp b/SU2_CFD/src/solver_direct_turbulent.cpp index 33d310d7dbd7..a92e8b8f03c5 100644 --- a/SU2_CFD/src/solver_direct_turbulent.cpp +++ b/SU2_CFD/src/solver_direct_turbulent.cpp @@ -971,8 +971,7 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (SA model)." << endl; Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } @@ -3311,8 +3310,7 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (SST model)." << endl; Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); - if ((config->GetKind_Linear_Solver_Prec() == LINELET) || - (config->GetKind_Linear_Solver() == SMOOTHER_LINELET)) { + if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); if (rank == MASTER_NODE) cout << "Compute linelet structure. " << nLineLets << " elements in each line (average)." << endl; } diff --git a/TestCases/aeroelastic/aeroelastic_NACA64A010.cfg b/TestCases/aeroelastic/aeroelastic_NACA64A010.cfg index 80acc8f90bc1..0e9763908888 100644 --- a/TestCases/aeroelastic/aeroelastic_NACA64A010.cfg +++ b/TestCases/aeroelastic/aeroelastic_NACA64A010.cfg @@ -171,9 +171,7 @@ CFL_NUMBER= 4.0 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/cont_adj_euler/naca0012/inv_NACA0012.cfg b/TestCases/cont_adj_euler/naca0012/inv_NACA0012.cfg index 2721106bf3bc..757d7f856a5f 100644 --- a/TestCases/cont_adj_euler/naca0012/inv_NACA0012.cfg +++ b/TestCases/cont_adj_euler/naca0012/inv_NACA0012.cfg @@ -106,9 +106,7 @@ EXT_ITER= 150 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/cont_adj_euler/naca0012/inv_NACA0012_FD.cfg b/TestCases/cont_adj_euler/naca0012/inv_NACA0012_FD.cfg index d36933dfd50b..659c0e76aa5f 100644 --- a/TestCases/cont_adj_euler/naca0012/inv_NACA0012_FD.cfg +++ b/TestCases/cont_adj_euler/naca0012/inv_NACA0012_FD.cfg @@ -107,9 +107,7 @@ EXT_ITER= 10 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/deformation/cylindrical_ffd/def_cylindrical.cfg b/TestCases/deformation/cylindrical_ffd/def_cylindrical.cfg index 519a574a0e52..76a802b19b6c 100644 --- a/TestCases/deformation/cylindrical_ffd/def_cylindrical.cfg +++ b/TestCases/deformation/cylindrical_ffd/def_cylindrical.cfg @@ -123,9 +123,7 @@ LIMITER_ITER= 999999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/deformation/naca4412/def_NACA4412.cfg b/TestCases/deformation/naca4412/def_NACA4412.cfg index 3abd0d9d8a86..bdab44e1cdaa 100644 --- a/TestCases/deformation/naca4412/def_NACA4412.cfg +++ b/TestCases/deformation/naca4412/def_NACA4412.cfg @@ -115,9 +115,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/deformation/spherical_ffd/def_spherical.cfg b/TestCases/deformation/spherical_ffd/def_spherical.cfg index c041eff5ae51..858ad51ca39d 100644 --- a/TestCases/deformation/spherical_ffd/def_spherical.cfg +++ b/TestCases/deformation/spherical_ffd/def_spherical.cfg @@ -123,9 +123,7 @@ LIMITER_ITER= 999999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/deformation/spherical_ffd/def_spherical_bspline.cfg b/TestCases/deformation/spherical_ffd/def_spherical_bspline.cfg index 361cb5f80100..83960aba71d3 100644 --- a/TestCases/deformation/spherical_ffd/def_spherical_bspline.cfg +++ b/TestCases/deformation/spherical_ffd/def_spherical_bspline.cfg @@ -123,9 +123,7 @@ LIMITER_ITER= 999999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg b/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg index e5ea67755481..75a1d7222629 100644 --- a/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg +++ b/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg @@ -360,9 +360,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/disc_adj_euler/cylinder3D/inv_cylinder3D.cfg b/TestCases/disc_adj_euler/cylinder3D/inv_cylinder3D.cfg index 77c4efa35c90..4d11eb702860 100644 --- a/TestCases/disc_adj_euler/cylinder3D/inv_cylinder3D.cfg +++ b/TestCases/disc_adj_euler/cylinder3D/inv_cylinder3D.cfg @@ -205,9 +205,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/disc_adj_fsi/Airfoil_2d/grad_young.opt.ref b/TestCases/disc_adj_fsi/Airfoil_2d/grad_young.opt.ref new file mode 100644 index 000000000000..531c86a9fe38 --- /dev/null +++ b/TestCases/disc_adj_fsi/Airfoil_2d/grad_young.opt.ref @@ -0,0 +1,2 @@ +INDEX GRAD +0 -1.712659639552295e-03 diff --git a/TestCases/disc_adj_incomp_navierstokes/cylinder/heated_cylinder.cfg b/TestCases/disc_adj_incomp_navierstokes/cylinder/heated_cylinder.cfg index 6ed4d856c9d1..6eaf44d6eccb 100644 --- a/TestCases/disc_adj_incomp_navierstokes/cylinder/heated_cylinder.cfg +++ b/TestCases/disc_adj_incomp_navierstokes/cylinder/heated_cylinder.cfg @@ -212,9 +212,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sa.cfg b/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sa.cfg index 980ad5aae39d..fd8d45f207bc 100644 --- a/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sa.cfg +++ b/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sa.cfg @@ -118,9 +118,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sst.cfg b/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sst.cfg index e24476ce7c6c..5b0a31c28c8b 100644 --- a/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sst.cfg +++ b/TestCases/disc_adj_rans/naca0012/turb_NACA0012_sst.cfg @@ -118,9 +118,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/disc_adj_turbomachinery/transonic_stator_2D/transonic_stator.cfg b/TestCases/disc_adj_turbomachinery/transonic_stator_2D/transonic_stator.cfg index f40fa1a3d967..903818991e22 100644 --- a/TestCases/disc_adj_turbomachinery/transonic_stator_2D/transonic_stator.cfg +++ b/TestCases/disc_adj_turbomachinery/transonic_stator_2D/transonic_stator.cfg @@ -212,7 +212,7 @@ CFL_ADAPT_PARAM= ( 1.3, 1.2, 1.0, 10.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/fea_topology/grad_ref_node.dat.ref b/TestCases/fea_topology/grad_ref_node.dat.ref index 6f3a4d4a6854..15afff813c19 100644 --- a/TestCases/fea_topology/grad_ref_node.dat.ref +++ b/TestCases/fea_topology/grad_ref_node.dat.ref @@ -6157,7 +6157,7 @@ -9.57383e-11 -2.60763e-10 -2.54078e-10 --1.57598e-10 +-1.57599e-10 -8.68931e-12 -4.26927e-13 0 diff --git a/TestCases/gust/inv_gust_NACA0012.cfg b/TestCases/gust/inv_gust_NACA0012.cfg index a1e781476de8..8ffe172b5c85 100644 --- a/TestCases/gust/inv_gust_NACA0012.cfg +++ b/TestCases/gust/inv_gust_NACA0012.cfg @@ -160,9 +160,7 @@ SENS_REMOVE_SHARP= NO % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/incomp_euler/nozzle/inv_nozzle.cfg b/TestCases/incomp_euler/nozzle/inv_nozzle.cfg index b42fa33f2e2a..1c7a7c5d1dce 100644 --- a/TestCases/incomp_euler/nozzle/inv_nozzle.cfg +++ b/TestCases/incomp_euler/nozzle/inv_nozzle.cfg @@ -115,9 +115,7 @@ CFL_NUMBER= 15.0 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/incomp_navierstokes/cylinder/poly_cylinder.cfg b/TestCases/incomp_navierstokes/cylinder/poly_cylinder.cfg index a771f19e710e..49338dcbf969 100644 --- a/TestCases/incomp_navierstokes/cylinder/poly_cylinder.cfg +++ b/TestCases/incomp_navierstokes/cylinder/poly_cylinder.cfg @@ -193,9 +193,7 @@ JST_SENSOR_COEFF= ( 0.5, 0.02 ) % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index f8571bb56b74..0c492033e785 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -866,7 +866,7 @@ def main(): channel_2D.cfg_dir = "sliding_interface/channel_2D" channel_2D.cfg_file = "channel_2D_WA.cfg" channel_2D.test_iter = 4 - channel_2D.test_vals = [-1.656863, 4.263134, 0.000000, 0.000000] #last 4 columns + channel_2D.test_vals = [-1.656871, 4.263104, 0.000000, 0.000000] #last 4 columns channel_2D.su2_exec = "parallel_computation.py -f" channel_2D.timeout = 100 channel_2D.tol = 0.00001 @@ -902,7 +902,7 @@ def main(): rotating_cylinders.cfg_dir = "sliding_interface/rotating_cylinders" rotating_cylinders.cfg_file = "rot_cylinders_WA.cfg" rotating_cylinders.test_iter = 3 - rotating_cylinders.test_vals = [1.219987, 7.729743, 0.000000, 0.000000] #last 4 columns + rotating_cylinders.test_vals = [-1.254581, 4.531726, 0.000000, 0.000000] #last 4 columns rotating_cylinders.su2_exec = "parallel_computation.py -f" rotating_cylinders.timeout = 1600 rotating_cylinders.tol = 0.00001 @@ -914,7 +914,7 @@ def main(): supersonic_vortex_shedding.cfg_dir = "sliding_interface/supersonic_vortex_shedding" supersonic_vortex_shedding.cfg_file = "sup_vor_shed_WA.cfg" supersonic_vortex_shedding.test_iter = 5 - supersonic_vortex_shedding.test_vals = [-1.124318, 4.605281, 0.000000, 0.000000] #last 4 columns + supersonic_vortex_shedding.test_vals = [-1.124330, 4.605307, 0.000000, 0.000000] #last 4 columns supersonic_vortex_shedding.su2_exec = "parallel_computation.py -f" supersonic_vortex_shedding.timeout = 1600 supersonic_vortex_shedding.tol = 0.00001 @@ -964,7 +964,7 @@ def main(): statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" statbeam3d.cfg_file = "configBeam_3d.cfg" statbeam3d.test_iter = 0 - statbeam3d.test_vals = [-8.396805, -8.162227, -8.156107, 64095.000000] #last 4 columns + statbeam3d.test_vals = [-8.396797, -8.162206, -8.156102, 64095.0] #last 4 columns statbeam3d.su2_exec = "parallel_computation_fsi.py -f" statbeam3d.timeout = 1600 statbeam3d.tol = 0.00001 diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 4be85b76ba62..bc5aecb9002f 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -223,7 +223,7 @@ def main(): discadj_fea.cfg_dir = "disc_adj_fea" discadj_fea.cfg_file = "configAD_fem.cfg" discadj_fea.test_iter = 9 - discadj_fea.test_vals = [-4.751576, -4.692029, -0.000364, -8.708700] #last 4 columns + discadj_fea.test_vals = [-5.394766, -5.572142, -0.000364, -8.708700] #last 4 columns discadj_fea.su2_exec = "parallel_computation.py -f" discadj_fea.timeout = 1600 discadj_fea.tol = 0.00001 diff --git a/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg b/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg index ba251e577267..b3f01df5408d 100644 --- a/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg +++ b/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg @@ -390,9 +390,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/py_wrapper/flatPlate_unsteady_CHT/unsteady_CHT_FlatPlate_Conf.cfg b/TestCases/py_wrapper/flatPlate_unsteady_CHT/unsteady_CHT_FlatPlate_Conf.cfg index 4d8051718181..353cd6c785a7 100644 --- a/TestCases/py_wrapper/flatPlate_unsteady_CHT/unsteady_CHT_FlatPlate_Conf.cfg +++ b/TestCases/py_wrapper/flatPlate_unsteady_CHT/unsteady_CHT_FlatPlate_Conf.cfg @@ -390,9 +390,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans/naca0012/turb_NACA0012_sa.cfg b/TestCases/rans/naca0012/turb_NACA0012_sa.cfg index 3bdabc984197..a889acc2f288 100644 --- a/TestCases/rans/naca0012/turb_NACA0012_sa.cfg +++ b/TestCases/rans/naca0012/turb_NACA0012_sa.cfg @@ -114,9 +114,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans/naca0012/turb_NACA0012_sa_binary.cfg b/TestCases/rans/naca0012/turb_NACA0012_sa_binary.cfg index ed1df3d2450e..08479a405849 100644 --- a/TestCases/rans/naca0012/turb_NACA0012_sa_binary.cfg +++ b/TestCases/rans/naca0012/turb_NACA0012_sa_binary.cfg @@ -120,9 +120,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans/naca0012/turb_NACA0012_sa_neg.cfg b/TestCases/rans/naca0012/turb_NACA0012_sa_neg.cfg index 2c20dd607ddc..0413ab0545fc 100644 --- a/TestCases/rans/naca0012/turb_NACA0012_sa_neg.cfg +++ b/TestCases/rans/naca0012/turb_NACA0012_sa_neg.cfg @@ -108,9 +108,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans/naca0012/turb_NACA0012_sst.cfg b/TestCases/rans/naca0012/turb_NACA0012_sst.cfg index b4dc82191d6e..f532c0c5f6fc 100644 --- a/TestCases/rans/naca0012/turb_NACA0012_sst.cfg +++ b/TestCases/rans/naca0012/turb_NACA0012_sst.cfg @@ -114,9 +114,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans/naca0012/turb_NACA0012_sst_multigrid_restart.cfg b/TestCases/rans/naca0012/turb_NACA0012_sst_multigrid_restart.cfg index 8f10b13f391b..40cb1d447976 100644 --- a/TestCases/rans/naca0012/turb_NACA0012_sst_multigrid_restart.cfg +++ b/TestCases/rans/naca0012/turb_NACA0012_sst_multigrid_restart.cfg @@ -247,9 +247,7 @@ SENS_REMOVE_SHARP= NO % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans/propeller/propeller.cfg b/TestCases/rans/propeller/propeller.cfg index 4bfec0fe08a1..1f72b1142610 100644 --- a/TestCases/rans/propeller/propeller.cfg +++ b/TestCases/rans/propeller/propeller.cfg @@ -151,9 +151,7 @@ OBJECTIVE_FUNCTION= DRAG % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans_uq/naca0012/turb_NACA0012_uq.cfg b/TestCases/rans_uq/naca0012/turb_NACA0012_uq.cfg index d8b06d80efa1..a2faaa2b3541 100644 --- a/TestCases/rans_uq/naca0012/turb_NACA0012_uq.cfg +++ b/TestCases/rans_uq/naca0012/turb_NACA0012_uq.cfg @@ -105,9 +105,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_1c.cfg b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_1c.cfg index 8fd6df48f453..2c908bb33e86 100644 --- a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_1c.cfg +++ b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_1c.cfg @@ -122,9 +122,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_2c.cfg b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_2c.cfg index 84a1a83c9030..7ef9cf87a4b8 100644 --- a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_2c.cfg +++ b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_2c.cfg @@ -122,9 +122,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_3c.cfg b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_3c.cfg index e3fe59840638..eb4f91e6fa08 100644 --- a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_3c.cfg +++ b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_3c.cfg @@ -122,9 +122,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_p1c1.cfg b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_p1c1.cfg index 60592902876e..9fb184fce76c 100644 --- a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_p1c1.cfg +++ b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_p1c1.cfg @@ -122,9 +122,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_p1c2.cfg b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_p1c2.cfg index f5f6156fd744..0cddf04b725d 100644 --- a/TestCases/rans_uq/naca0012/turb_NACA0012_uq_p1c2.cfg +++ b/TestCases/rans_uq/naca0012/turb_NACA0012_uq_p1c2.cfg @@ -122,9 +122,7 @@ LIMITER_ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, MULTIGRID, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 3b04fe4c9674..00987391d19c 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -877,7 +877,7 @@ def main(): channel_2D.cfg_dir = "sliding_interface/channel_2D" channel_2D.cfg_file = "channel_2D_WA.cfg" channel_2D.test_iter = 4 - channel_2D.test_vals = [-1.656852, 4.263152, 0.000000, 0.000000] #last 4 columns + channel_2D.test_vals = [-1.656891, 4.263153, 0.000000, 0.000000] #last 4 columns channel_2D.su2_exec = "SU2_CFD" channel_2D.timeout = 100 channel_2D.tol = 0.00001 @@ -925,7 +925,7 @@ def main(): supersonic_vortex_shedding.cfg_dir = "sliding_interface/supersonic_vortex_shedding" supersonic_vortex_shedding.cfg_file = "sup_vor_shed_WA.cfg" supersonic_vortex_shedding.test_iter = 5 - supersonic_vortex_shedding.test_vals = [-1.130591, 4.595041, 0.000000, 0.000000] #last 4 columns + supersonic_vortex_shedding.test_vals = [-1.129509, 4.596356, 0.000000, 0.000000] #last 4 columns supersonic_vortex_shedding.su2_exec = "SU2_CFD" supersonic_vortex_shedding.timeout = 1600 supersonic_vortex_shedding.tol = 0.00001 @@ -975,7 +975,7 @@ def main(): statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" statbeam3d.cfg_file = "configBeam_3d.cfg" statbeam3d.test_iter = 0 - statbeam3d.test_vals = [-8.498274, -8.230638, -8.123824, 6.4095e+04] #last 4 columns + statbeam3d.test_vals = [-8.498245, -8.230816, -8.123810, 64095.0] #last 4 columns statbeam3d.su2_exec = "SU2_CFD" statbeam3d.timeout = 1600 statbeam3d.tol = 0.00001 @@ -1019,7 +1019,7 @@ def main(): airfoilRBF.cfg_dir = "fea_fsi/Airfoil_RBF" airfoilRBF.cfg_file = "config.cfg" airfoilRBF.test_iter = 29 - airfoilRBF.test_vals = [-13.474206, -3.771392, -11.646689, 1.4045e+06] #last 4 columns + airfoilRBF.test_vals = [-13.228126, -3.738614, -11.639760, 1404500.0] #last 4 columns airfoilRBF.su2_exec = "SU2_CFD" airfoilRBF.timeout = 1600 airfoilRBF.tol = 0.00001 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 394518c89357..90dda349e6ee 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -208,7 +208,7 @@ def main(): discadj_fea.cfg_dir = "disc_adj_fea" discadj_fea.cfg_file = "configAD_fem.cfg" discadj_fea.test_iter = 9 - discadj_fea.test_vals = [-6.282767, -6.361594, -3.6413e-04, -8.7087e+00] #last 4 columns + discadj_fea.test_vals = [-6.319841, -6.375512, -0.000364, -8.708700] #last 4 columns discadj_fea.su2_exec = "SU2_CFD_AD" discadj_fea.timeout = 1600 discadj_fea.tol = 0.00001 diff --git a/TestCases/sliding_interface/channel_2D/channel_2D_NN.cfg b/TestCases/sliding_interface/channel_2D/channel_2D_NN.cfg index 3183ca68602c..26a5376b353b 100644 --- a/TestCases/sliding_interface/channel_2D/channel_2D_NN.cfg +++ b/TestCases/sliding_interface/channel_2D/channel_2D_NN.cfg @@ -122,7 +122,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/channel_2D/channel_2D_WA.cfg b/TestCases/sliding_interface/channel_2D/channel_2D_WA.cfg index b741b61fd2d8..53a0488edcaf 100644 --- a/TestCases/sliding_interface/channel_2D/channel_2D_WA.cfg +++ b/TestCases/sliding_interface/channel_2D/channel_2D_WA.cfg @@ -123,7 +123,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/channel_3D/channel_3D_NN.cfg b/TestCases/sliding_interface/channel_3D/channel_3D_NN.cfg index 09b3334514a5..8cbfcb5d7a1a 100644 --- a/TestCases/sliding_interface/channel_3D/channel_3D_NN.cfg +++ b/TestCases/sliding_interface/channel_3D/channel_3D_NN.cfg @@ -124,7 +124,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/channel_3D/channel_3D_WA.cfg b/TestCases/sliding_interface/channel_3D/channel_3D_WA.cfg index cd0feb0ede6f..b092b9f002d0 100644 --- a/TestCases/sliding_interface/channel_3D/channel_3D_WA.cfg +++ b/TestCases/sliding_interface/channel_3D/channel_3D_WA.cfg @@ -125,7 +125,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/pipe/pipe_NN.cfg b/TestCases/sliding_interface/pipe/pipe_NN.cfg index f5254531ac3e..52c5a1b82821 100644 --- a/TestCases/sliding_interface/pipe/pipe_NN.cfg +++ b/TestCases/sliding_interface/pipe/pipe_NN.cfg @@ -127,7 +127,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/pipe/pipe_WA.cfg b/TestCases/sliding_interface/pipe/pipe_WA.cfg index 09743629389e..6572068cb2b5 100644 --- a/TestCases/sliding_interface/pipe/pipe_WA.cfg +++ b/TestCases/sliding_interface/pipe/pipe_WA.cfg @@ -128,7 +128,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/rotating_cylinders/rot_cylinders_NN.cfg b/TestCases/sliding_interface/rotating_cylinders/rot_cylinders_NN.cfg index 49eb02a8df05..bb8f6e2d769b 100644 --- a/TestCases/sliding_interface/rotating_cylinders/rot_cylinders_NN.cfg +++ b/TestCases/sliding_interface/rotating_cylinders/rot_cylinders_NN.cfg @@ -128,7 +128,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/rotating_cylinders/rot_cylinders_WA.cfg b/TestCases/sliding_interface/rotating_cylinders/rot_cylinders_WA.cfg index 6f8092bad311..4e9811441731 100644 --- a/TestCases/sliding_interface/rotating_cylinders/rot_cylinders_WA.cfg +++ b/TestCases/sliding_interface/rotating_cylinders/rot_cylinders_WA.cfg @@ -129,7 +129,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/single_stage/single_stage_NN.cfg b/TestCases/sliding_interface/single_stage/single_stage_NN.cfg index 6bab2bf9e5b5..ddf2c3864fd8 100644 --- a/TestCases/sliding_interface/single_stage/single_stage_NN.cfg +++ b/TestCases/sliding_interface/single_stage/single_stage_NN.cfg @@ -142,7 +142,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/single_stage/single_stage_WA.cfg b/TestCases/sliding_interface/single_stage/single_stage_WA.cfg index bbfeca961c5e..f8df591d109a 100644 --- a/TestCases/sliding_interface/single_stage/single_stage_WA.cfg +++ b/TestCases/sliding_interface/single_stage/single_stage_WA.cfg @@ -142,7 +142,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/supersonic_vortex_shedding/sup_vor_shed_NN.cfg b/TestCases/sliding_interface/supersonic_vortex_shedding/sup_vor_shed_NN.cfg index 14b6372b3bcb..63851b2170f6 100644 --- a/TestCases/sliding_interface/supersonic_vortex_shedding/sup_vor_shed_NN.cfg +++ b/TestCases/sliding_interface/supersonic_vortex_shedding/sup_vor_shed_NN.cfg @@ -129,7 +129,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/supersonic_vortex_shedding/sup_vor_shed_WA.cfg b/TestCases/sliding_interface/supersonic_vortex_shedding/sup_vor_shed_WA.cfg index f7ce85894625..b978a254bcf5 100644 --- a/TestCases/sliding_interface/supersonic_vortex_shedding/sup_vor_shed_WA.cfg +++ b/TestCases/sliding_interface/supersonic_vortex_shedding/sup_vor_shed_WA.cfg @@ -130,7 +130,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/uniform_flow/uniform_NN.cfg b/TestCases/sliding_interface/uniform_flow/uniform_NN.cfg index a38f8e73d0c6..172b67d5e424 100644 --- a/TestCases/sliding_interface/uniform_flow/uniform_NN.cfg +++ b/TestCases/sliding_interface/uniform_flow/uniform_NN.cfg @@ -129,7 +129,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/sliding_interface/uniform_flow/uniform_WA.cfg b/TestCases/sliding_interface/uniform_flow/uniform_WA.cfg index abc3adcb4202..4f2504686d1c 100644 --- a/TestCases/sliding_interface/uniform_flow/uniform_WA.cfg +++ b/TestCases/sliding_interface/uniform_flow/uniform_WA.cfg @@ -130,7 +130,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/turbomachinery/APU_turbocharger/Jones.cfg b/TestCases/turbomachinery/APU_turbocharger/Jones.cfg index f98025a5b4e5..2d882b487afc 100755 --- a/TestCases/turbomachinery/APU_turbocharger/Jones.cfg +++ b/TestCases/turbomachinery/APU_turbocharger/Jones.cfg @@ -213,7 +213,7 @@ CFL_ADAPT_PARAM= ( 1.3, 1.2, 5.0, 30.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/turbomachinery/APU_turbocharger/Jones_rst.cfg b/TestCases/turbomachinery/APU_turbocharger/Jones_rst.cfg index 5528a7a2c774..f9f075dc6ad2 100755 --- a/TestCases/turbomachinery/APU_turbocharger/Jones_rst.cfg +++ b/TestCases/turbomachinery/APU_turbocharger/Jones_rst.cfg @@ -218,7 +218,7 @@ CFL_ADAPT_PARAM= ( 1.3, 1.2, 5.0, 30.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/turbomachinery/axial_stage_2D/Axial_stage2D.cfg b/TestCases/turbomachinery/axial_stage_2D/Axial_stage2D.cfg index bcb864b7a7b9..091091220e91 100755 --- a/TestCases/turbomachinery/axial_stage_2D/Axial_stage2D.cfg +++ b/TestCases/turbomachinery/axial_stage_2D/Axial_stage2D.cfg @@ -203,7 +203,7 @@ TRANSLATION_RATE_Y= 0.0 -150.0 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/turbomachinery/centrifugal_blade/centrifugal_blade.cfg b/TestCases/turbomachinery/centrifugal_blade/centrifugal_blade.cfg index 0c5066ae4fca..ce3133fc6253 100755 --- a/TestCases/turbomachinery/centrifugal_blade/centrifugal_blade.cfg +++ b/TestCases/turbomachinery/centrifugal_blade/centrifugal_blade.cfg @@ -205,7 +205,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.4, 10.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/turbomachinery/centrifugal_stage/centrifugal_stage.cfg b/TestCases/turbomachinery/centrifugal_stage/centrifugal_stage.cfg index 247e000d4b08..179ad4bc2c99 100755 --- a/TestCases/turbomachinery/centrifugal_stage/centrifugal_stage.cfg +++ b/TestCases/turbomachinery/centrifugal_stage/centrifugal_stage.cfg @@ -212,7 +212,7 @@ CFL_ADAPT_PARAM= ( 0.3, 0.5, 1.0, 1000.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/turbomachinery/transonic_stator_2D/transonic_stator.cfg b/TestCases/turbomachinery/transonic_stator_2D/transonic_stator.cfg index 06c2f86c0347..75e7b355548d 100644 --- a/TestCases/turbomachinery/transonic_stator_2D/transonic_stator.cfg +++ b/TestCases/turbomachinery/transonic_stator_2D/transonic_stator.cfg @@ -201,7 +201,7 @@ CFL_ADAPT_PARAM= ( 1.3, 1.2, 1.0, 10.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/TestCases/turbomachinery/transonic_stator_2D/transonic_stator_rst.cfg b/TestCases/turbomachinery/transonic_stator_2D/transonic_stator_rst.cfg index e313fe53f200..3ac2beec2064 100644 --- a/TestCases/turbomachinery/transonic_stator_2D/transonic_stator_rst.cfg +++ b/TestCases/turbomachinery/transonic_stator_2D/transonic_stator_rst.cfg @@ -206,7 +206,7 @@ CFL_ADAPT_PARAM= ( 1.3, 1.2, 1.0, 10.0) % % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, SMOOTHER_ILU, SMOOTHER_LUSGS, SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES % % Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) diff --git a/config_template.cfg b/config_template.cfg index 017538116113..680cfea489f1 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -12,7 +12,7 @@ % ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% % % Physical governing equations (EULER, NAVIER_STOKES, - FEM_EULER, FEM_NAVIER_STOKES, FEM_RANS, FEM_LES, +% FEM_EULER, FEM_NAVIER_STOKES, FEM_RANS, FEM_LES, % WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, % POISSON_EQUATION) PHYSICAL_PROBLEM= EULER @@ -843,14 +843,19 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER_JACOBI, -% SMOOTHER_ILU, SMOOTHER_LUSGS, -% SMOOTHER_LINELET) +% Linear solver or smoother for implicit formulations: +% BCGSTAB, FGMRES, RESTARTED_FGMRES, CONJUGATE_GRADIENT (self-adjoint problems only), SMOOTHER. LINEAR_SOLVER= FGMRES % -% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) +% Same for discrete adjoint (smoothers not supported) +DISCADJ_LIN_SOLVER= FGMRES +% +% Preconditioner of the Krylov linear solver or type of smoother (ILU, LU_SGS, LINELET, JACOBI) LINEAR_SOLVER_PREC= ILU % +% Same for discrete adjoint (JACOBI or ILU) +DISCADJ_LIN_PREC= ILU +% % Linael solver ILU preconditioner fill-in level (0 by default) LINEAR_SOLVER_ILU_FILL_IN= 0 % @@ -859,6 +864,12 @@ LINEAR_SOLVER_ERROR= 1E-6 % % Max number of iterations of the linear solver for the implicit formulation LINEAR_SOLVER_ITER= 5 +% +% Restart frequency for RESTARTED_FGMRES +LINEAR_SOLVER_RESTART_FREQUENCY= 10 +% +% Relaxation factor for smoother-type solvers (LINEAR_SOLVER= SMOOTHER) +LINEAR_SOLVER_SMOOTHER_RELAXATION= 1.0 % -------------------------- MULTIGRID PARAMETERS -----------------------------% %