diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index 82a3b92b101c..ecd9e069edd0 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -16,7 +16,7 @@ jobs: strategy: fail-fast: false matrix: - config_set: [BaseMPI, ReverseMPI, ForwardMPI, BaseNoMPI, ReverseNoMPI, ForwardNoMPI] + config_set: [BaseMPI, ReverseMPI, ForwardMPI, BaseNoMPI, ReverseNoMPI, ForwardNoMPI, BaseOMP] include: - config_set: BaseMPI flags: '-Denable-pywrapper=true -Denable-tests=true --werror' @@ -30,6 +30,8 @@ jobs: flags: '-Denable-autodiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-pywrapper=true -Denable-tests=true --werror' - config_set: ForwardNoMPI flags: '-Denable-directdiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-tests=true --werror' + - config_set: BaseOMP + flags: '-Dwith-omp=true -Denable-mixedprec=true -Denable-tecio=false --werror' runs-on: ubuntu-latest steps: - name: Cache Object Files @@ -54,7 +56,7 @@ jobs: strategy: fail-fast: false matrix: - testscript: ['tutorials.py', 'parallel_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py'] + testscript: ['tutorials.py', 'parallel_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py', 'hybrid_regression.py'] include: - testscript: 'tutorials.py' tag: MPI @@ -66,31 +68,22 @@ jobs: tag: NoMPI - testscript: 'serial_regression_AD.py' tag: NoMPI + - testscript: 'hybrid_regression.py' + tag: OMP steps: - - name: Download Base - uses: actions/download-artifact@v1 - with: - name: ${{format('Base{0}', matrix.tag)}} - - name: Download Reverse - uses: actions/download-artifact@v1 - with: - name: ${{format('Reverse{0}', matrix.tag)}} - - name: Download Forward - uses: actions/download-artifact@v1 - with: - name: ${{format('Forward{0}', matrix.tag)}} + - name: Download All artifact + uses: actions/download-artifact@v2 - name: Move Binaries run: | mkdir -p install/bin - cp -r ${{format('Base{0}', matrix.tag)}}/. install/bin/ - cp -r ${{format('Reverse{0}', matrix.tag)}}/. install/bin/ - cp -r ${{format('Forward{0}', matrix.tag)}}/. install/bin/ - chmod a+x install/bin/* + if [ -d "${{format('Base{0}', matrix.tag)}}" ]; then cp -r ${{format('Base{0}', matrix.tag)}}/. install/bin/; fi + if [ -d "${{format('Reverse{0}', matrix.tag)}}" ]; then cp -r ${{format('Reverse{0}', matrix.tag)}}/. install/bin/; fi + if [ -d "${{format('Forward{0}', matrix.tag)}}" ]; then cp -r ${{format('Forward{0}', matrix.tag)}}/. install/bin/; fi + chmod a+x install/bin/* - name: Run Tests in Container uses: docker://su2code/test-su2:20200303 with: args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} - unit_tests: runs-on: ubuntu-latest name: Unit Tests @@ -125,7 +118,7 @@ jobs: cp -r ${{format('Base{0}', matrix.tag)}}/. install/bin/ cp -r ${{format('Reverse{0}', matrix.tag)}}/. install/bin/ cp -r ${{format('Forward{0}', matrix.tag)}}/. install/bin/ - chmod a+x install/bin/* + chmod a+x install/bin/* - name: Run Unit Tests uses: docker://su2code/test-su2:20191031 with: diff --git a/.travis.yml b/.travis.yml index 99dabda9cc4c..284818490731 100644 --- a/.travis.yml +++ b/.travis.yml @@ -37,11 +37,15 @@ env: - CONFIGURE_COMMAND="./meson.py build --prefix=$TRAVIS_BUILD_DIR -Denable-pywrapper=true" TEST_SCRIPT=parallel_regression.py + # Hybrid-parallel build and test + - CONFIGURE_COMMAND="./meson.py build --prefix=$TRAVIS_BUILD_DIR -Dwith-omp=true -Denable-mixedprec=true -Denable-tecio=false" + TEST_SCRIPT=hybrid_regression.py + # Serial build and test for AD - CONFIGURE_COMMAND="./meson.py build --prefix=$TRAVIS_BUILD_DIR -Denable-pywrapper=true -Dwith-mpi=disabled -Denable-autodiff=true -Denable-directdiff=true" TEST_SCRIPT=serial_regression_AD.py - # Parallel build and test for AD: + # Parallel build and test for AD - CONFIGURE_COMMAND="./meson.py build --prefix=$TRAVIS_BUILD_DIR -Denable-pywrapper=true -Denable-autodiff=true" TEST_SCRIPT=parallel_regression_AD.py @@ -51,7 +55,7 @@ before_install: - sudo update-alternatives --install /usr/bin/python python /usr/bin/python3 10 # Install the necessary packages using apt-get with sudo - sudo apt-get update -qq - - sudo apt-get install -qq build-essential python3-numpy python3-scipy libopenmpi-dev openmpi-bin swig python3-mpi4py + - sudo apt-get install -qq build-essential python3-numpy python3-scipy libopenmpi-dev openmpi-bin swig python3-mpi4py # to avoid interference with MPI - test -n $CC && unset CC - test -n $CXX && unset CXX diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index fbc668200cf0..3c7f70d33557 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -44,7 +44,6 @@ #include #include "./option_structure.hpp" -#include "./datatype_structure.hpp" #ifdef HAVE_CGNS #include "cgnslib.h" @@ -1282,7 +1281,7 @@ class CConfig { * \brief Constructor of the class which takes an istream buffer containing the config options. */ CConfig(istream &case_buffer, unsigned short val_software, bool verb_high); - + /*! * \brief Constructor of the class which reads the input file and uses default options from another config. */ @@ -1302,7 +1301,7 @@ class CConfig { * \brief Destructor of the class. */ ~CConfig(void); - + /*! * \brief Initialize common fields of the config structure. */ @@ -8009,8 +8008,8 @@ class CConfig { /*! * \brief Set the config file parsing. */ - void SetConfig_Parsing(istream &config_buffer); - + void SetConfig_Parsing(istream &config_buffer); + /*! * \brief Set the config file parsing. */ @@ -8340,26 +8339,6 @@ class CConfig { */ unsigned long GetNonphysical_Reconstr(void) const { return Nonphys_Reconstr; } - /*! - * \brief Given arrays x[1..n] and y[1..n] containing a tabulated function, i.e., yi = f(xi), with - x1 < x2 < . . . < xN , and given values yp1 and ypn for the first derivative of the interpolating - function at points 1 and n, respectively, this routine returns an array y2[1..n] that contains - the second derivatives of the interpolating function at the tabulated points xi. If yp1 and/or - ypn are equal to 1 × 1030 or larger, the routine is signaled to set the corresponding boundary - condition for a natural spline, with zero second derivative on that boundary. - Numerical Recipes: The Art of Scientific Computing, Third Edition in C++. - */ - void SetSpline(vector &x, vector &y, unsigned long n, su2double yp1, su2double ypn, vector &y2); - - /*! - * \brief Given the arrays xa[1..n] and ya[1..n], which tabulate a function (with the xai’s in order), - and given the array y2a[1..n], which is the output from spline above, and given a value of - x, this routine returns a cubic-spline interpolated value y. - Numerical Recipes: The Art of Scientific Computing, Third Edition in C++. - * \return The interpolated value of for x. - */ - su2double GetSpline(vector &xa, vector &ya, vector &y2a, unsigned long n, su2double x); - /*! * \brief Start the timer for profiling subroutines. * \param[in] val_start_time - the value of the start time. diff --git a/Common/include/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp similarity index 100% rename from Common/include/ad_structure.hpp rename to Common/include/basic_types/ad_structure.hpp diff --git a/Common/include/datatype_structure.hpp b/Common/include/basic_types/datatype_structure.hpp similarity index 100% rename from Common/include/datatype_structure.hpp rename to Common/include/basic_types/datatype_structure.hpp diff --git a/Common/include/blas_structure.hpp b/Common/include/blas_structure.hpp index 0d7c44673b0c..6cef28f83482 100644 --- a/Common/include/blas_structure.hpp +++ b/Common/include/blas_structure.hpp @@ -8,7 +8,7 @@ * * SU2 Project Website: https://su2code.github.io * - * The SU2 Project is maintained by the SU2 Foundation + * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md) @@ -29,7 +29,6 @@ #pragma once -#include "datatype_structure.hpp" #include "CConfig.hpp" /* LIBXSMM include files, if supported. */ @@ -39,7 +38,7 @@ /*! * \class CBlasStructure - * \brief Class, which serves as an interface to the BLAS functionalities needed. + * \brief Class, which serves as an interface to the BLAS functionalities needed. * \author: E. van der Weide * \version 7.0.4 "Blackbird" */ diff --git a/Common/include/fem/geometry_structure_fem_part.hpp b/Common/include/fem/geometry_structure_fem_part.hpp index fae54f7c009d..9b45bff487b3 100644 --- a/Common/include/fem/geometry_structure_fem_part.hpp +++ b/Common/include/fem/geometry_structure_fem_part.hpp @@ -27,7 +27,7 @@ #pragma once -#include "../datatype_structure.hpp" +#include "../basic_types/datatype_structure.hpp" #include #include diff --git a/Common/include/interface_interpolation/CInterpolator.hpp b/Common/include/interface_interpolation/CInterpolator.hpp index 574904265667..55cc6ba64331 100644 --- a/Common/include/interface_interpolation/CInterpolator.hpp +++ b/Common/include/interface_interpolation/CInterpolator.hpp @@ -26,7 +26,7 @@ */ #pragma once -#include "../../include/datatype_structure.hpp" +#include "../../include/basic_types/datatype_structure.hpp" #include "../../include/toolboxes/C2DContainer.hpp" #include diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 795a60802c59..c67117665f1a 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -353,7 +353,7 @@ class CSysMatrix { void Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar, unsigned short neqn, bool EdgeConnect, CGeometry *geometry, - CConfig *config, bool needTranspPtr = false); + const CConfig *config, bool needTranspPtr = false); /*! * \brief Sets to zero all the entries of the sparse matrix. diff --git a/Common/include/mpi_structure.hpp b/Common/include/mpi_structure.hpp index e8a024ba98a6..cb17027ad158 100644 --- a/Common/include/mpi_structure.hpp +++ b/Common/include/mpi_structure.hpp @@ -33,7 +33,7 @@ #include #endif -#include "./datatype_structure.hpp" +#include "basic_types/datatype_structure.hpp" #include #ifndef _MSC_VER #include diff --git a/Common/include/toolboxes/C1DInterpolation.hpp b/Common/include/toolboxes/C1DInterpolation.hpp index ede5236ef03a..1129a62c8f55 100644 --- a/Common/include/toolboxes/C1DInterpolation.hpp +++ b/Common/include/toolboxes/C1DInterpolation.hpp @@ -31,7 +31,6 @@ #include #include #include -#include "../datatype_structure.hpp" #include "../option_structure.hpp" using namespace std; diff --git a/Common/include/toolboxes/C2DContainer.hpp b/Common/include/toolboxes/C2DContainer.hpp index dcb8aac91392..cbc45f221d61 100644 --- a/Common/include/toolboxes/C2DContainer.hpp +++ b/Common/include/toolboxes/C2DContainer.hpp @@ -28,7 +28,7 @@ #pragma once #include "allocation_toolbox.hpp" -#include "../datatype_structure.hpp" +#include "../basic_types/datatype_structure.hpp" #include #include diff --git a/Common/include/toolboxes/printing_toolbox.hpp b/Common/include/toolboxes/printing_toolbox.hpp index dd1dd8ec07cc..c6280ae7a8be 100644 --- a/Common/include/toolboxes/printing_toolbox.hpp +++ b/Common/include/toolboxes/printing_toolbox.hpp @@ -32,7 +32,7 @@ #include #include #include -#include "../datatype_structure.hpp" +#include "../basic_types/datatype_structure.hpp" namespace PrintingToolbox { diff --git a/Common/lib/Makefile.am b/Common/lib/Makefile.am index 792bf2cf5351..e75362117462 100644 --- a/Common/lib/Makefile.am +++ b/Common/lib/Makefile.am @@ -59,7 +59,7 @@ lib_sources = \ ../src/grid_adaptation_structure.cpp \ ../src/grid_movement_structure.cpp \ ../src/mpi_structure.cpp \ - ../src/ad_structure.cpp \ + ../src/basic_types/ad_structure.cpp \ ../src/fem/fem_gauss_jacobi_quadrature.cpp \ ../src/geometry/CGeometry.cpp \ ../src/geometry/CPhysicalGeometry.cpp \ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 9ecbfd468f8c..f978ee8ba0aa 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -32,7 +32,7 @@ #include "../include/fem/fem_gauss_jacobi_quadrature.hpp" #include "../include/fem/fem_geometry_structure.hpp" -#include "../include/ad_structure.hpp" +#include "../include/basic_types/ad_structure.hpp" #include "../include/toolboxes/printing_toolbox.hpp" using namespace PrintingToolbox; @@ -3659,99 +3659,58 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ that for each option, a value has been declared for each moving marker. ---*/ if (nMarker_Moving > 0){ - unsigned short iDim; if (nMarkerMotion_Origin == 0){ nMarkerMotion_Origin = 3*nMarker_Moving; - MarkerMotion_Origin = new su2double[nMarkerMotion_Origin]; - for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ - for (iDim = 0; iDim < 3; iDim++){ - MarkerMotion_Origin[3*iMarker+iDim] = 0.0; - } - } + MarkerMotion_Origin = new su2double[nMarkerMotion_Origin] (); } if (nMarkerMotion_Origin/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_MOTION_ORIGIN must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); } if (nMarkerTranslation == 0){ nMarkerTranslation = 3*nMarker_Moving; - MarkerTranslation_Rate = new su2double[nMarkerTranslation]; - for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ - for (iDim = 0; iDim < 3; iDim++){ - MarkerTranslation_Rate[3*iMarker+iDim] = 0.0; - } - } + MarkerTranslation_Rate = new su2double[nMarkerTranslation] (); } if (nMarkerTranslation/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_TRANSLATION_RATE must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); } if (nMarkerRotation_Rate == 0){ nMarkerRotation_Rate = 3*nMarker_Moving; - MarkerRotation_Rate = new su2double[nMarkerRotation_Rate]; - for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ - for (iDim = 0; iDim < 3; iDim++){ - MarkerRotation_Rate[3*iMarker+iDim] = 0.0; - } - } + MarkerRotation_Rate = new su2double[nMarkerRotation_Rate] (); } if (nMarkerRotation_Rate/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_ROTATION_RATE must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); } if (nMarkerPlunging_Ampl == 0){ nMarkerPlunging_Ampl = 3*nMarker_Moving; - MarkerPlunging_Ampl = new su2double[nMarkerPlunging_Ampl]; - for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ - for (iDim = 0; iDim < 3; iDim++){ - MarkerPlunging_Ampl[3*iMarker+iDim] = 0.0; - } - } + MarkerPlunging_Ampl = new su2double[nMarkerPlunging_Ampl] (); } if (nMarkerPlunging_Ampl/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PLUNGING_AMPL must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); } if (nMarkerPlunging_Omega == 0){ nMarkerPlunging_Omega = 3*nMarker_Moving; - MarkerPlunging_Omega = new su2double[nMarkerPlunging_Omega]; - for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ - for (iDim = 0; iDim < 3; iDim++){ - MarkerPlunging_Omega[3*iMarker+iDim] = 0.0; - } - } + MarkerPlunging_Omega = new su2double[nMarkerPlunging_Omega] (); } if (nMarkerPlunging_Omega/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PLUNGING_OMEGA must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); } if (nMarkerPitching_Ampl == 0){ nMarkerPitching_Ampl = 3*nMarker_Moving; - MarkerPitching_Ampl = new su2double[nMarkerPitching_Ampl]; - for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ - for (iDim = 0; iDim < 3; iDim++){ - MarkerPitching_Ampl[3*iMarker+iDim] = 0.0; - } - } + MarkerPitching_Ampl = new su2double[nMarkerPitching_Ampl] (); } if (nMarkerPitching_Ampl/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PITCHING_AMPL must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); } if (nMarkerPitching_Omega == 0){ nMarkerPitching_Omega = 3*nMarker_Moving; - MarkerPitching_Omega = new su2double[nMarkerPitching_Omega]; - for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ - for (iDim = 0; iDim < 3; iDim++){ - MarkerPitching_Omega[3*iMarker+iDim] = 0.0; - } - } + MarkerPitching_Omega = new su2double[nMarkerPitching_Omega] (); } if (nMarkerPitching_Omega/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PITCHING_OMEGA must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); } if (nMarkerPitching_Phase == 0){ nMarkerPitching_Phase = 3*nMarker_Moving; - MarkerPitching_Phase = new su2double[nMarkerPitching_Phase]; - for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ - for (iDim = 0; iDim < 3; iDim++){ - MarkerPitching_Phase[3*iMarker+iDim] = 0.0; - } - } + MarkerPitching_Phase = new su2double[nMarkerPitching_Phase] (); } if (nMarkerPitching_Phase/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PITCHING_PHASE must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); @@ -3761,7 +3720,7 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ nMoveMotion_Origin = nMarker_Moving; MoveMotion_Origin = new unsigned short[nMoveMotion_Origin]; for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ - MoveMotion_Origin[iMarker] = NO; + MoveMotion_Origin[iMarker] = NO; } } if (nMoveMotion_Origin != nMarker_Moving){ @@ -9095,61 +9054,6 @@ short CConfig::FindInterfaceMarker(unsigned short iInterface) const { return -1; } -void CConfig::SetSpline(vector &x, vector &y, unsigned long n, su2double yp1, su2double ypn, vector &y2) { - unsigned long i, k; - su2double p, qn, sig, un, *u; - - u = new su2double [n]; - - if (yp1 > 0.99e30) // The lower boundary condition is set either to be "nat - y2[0]=u[0]=0.0; // -ural" - else { // or else to have a specified first derivative. - y2[0] = -0.5; - u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1); - } - - for (i=2; i<=n-1; i++) { // This is the decomposition loop of the tridiagonal al- - sig=(x[i-1]-x[i-2])/(x[i]-x[i-2]); // gorithm. y2 and u are used for tem- - p=sig*y2[i-2]+2.0; // porary storage of the decomposed - y2[i-1]=(sig-1.0)/p; // factors. - u[i-1]=(y[i]-y[i-1])/(x[i]-x[i-1]) - (y[i-1]-y[i-2])/(x[i-1]-x[i-2]); - u[i-1]=(6.0*u[i-1]/(x[i]-x[i-2])-sig*u[i-2])/p; - } - - if (ypn > 0.99e30) // The upper boundary condition is set either to be - qn=un=0.0; // "natural" - else { // or else to have a specified first derivative. - qn=0.5; - un=(3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2])); - } - y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0); - for (k=n-1; k>=1; k--) // This is the backsubstitution loop of the tridiagonal - y2[k-1]=y2[k-1]*y2[k]+u[k-1]; // algorithm. - - delete[] u; - -} - -su2double CConfig::GetSpline(vector&xa, vector&ya, vector&y2a, unsigned long n, su2double x) { - unsigned long klo, khi, k; - su2double h, b, a, y; - - klo=1; // We will find the right place in the table by means of - khi=n; // bisection. This is optimal if sequential calls to this - while (khi-klo > 1) { // routine are at random values of x. If sequential calls - k=(khi+klo) >> 1; // are in order, and closely spaced, one would do better - if (xa[k-1] > x) khi=k; // to store previous values of klo and khi and test if - else klo=k; // they remain appropriate on the next call. - } // klo and khi now bracket the input value of x - h=xa[khi-1]-xa[klo-1]; - if (h == 0.0) cout << "Bad xa input to routine splint" << endl; // The xa?s must be dis- - a=(xa[khi-1]-x)/h; // tinct. - b=(x-xa[klo-1])/h; // Cubic spline polynomial is now evaluated. - y=a*ya[klo-1]+b*ya[khi-1]+((a*a*a-a)*y2a[klo-1]+(b*b*b-b)*y2a[khi-1])*(h*h)/6.0; - - return y; -} - void CConfig::Tick(double *val_start_time) { #ifdef PROFILE diff --git a/Common/src/ad_structure.cpp b/Common/src/basic_types/ad_structure.cpp similarity index 96% rename from Common/src/ad_structure.cpp rename to Common/src/basic_types/ad_structure.cpp index bba75585d0e7..53d5a79e5104 100644 --- a/Common/src/ad_structure.cpp +++ b/Common/src/basic_types/ad_structure.cpp @@ -25,7 +25,7 @@ * License along with SU2. If not, see . */ -#include "../include/datatype_structure.hpp" +#include "../../include/basic_types/datatype_structure.hpp" namespace AD { #ifdef CODI_REVERSE_TYPE diff --git a/Common/src/fem/fem_cgns_elements.cpp b/Common/src/fem/fem_cgns_elements.cpp index fef883830506..ad0c7a8f2d8d 100644 --- a/Common/src/fem/fem_cgns_elements.cpp +++ b/Common/src/fem/fem_cgns_elements.cpp @@ -29,7 +29,6 @@ #include "../../include/fem/fem_cgns_elements.hpp" #include "../../include/fem/geometry_structure_fem_part.hpp" #include "../../include/mpi_structure.hpp" -#include "../../include/datatype_structure.hpp" #include #include diff --git a/Common/src/grid_movement_structure.cpp b/Common/src/grid_movement_structure.cpp index de9c64b7d92f..d280a58a8f11 100644 --- a/Common/src/grid_movement_structure.cpp +++ b/Common/src/grid_movement_structure.cpp @@ -195,7 +195,6 @@ void CVolumetricMovement::SetVolume_Deformation(CGeometry *geometry, CConfig *co } else if (Derivative && (config->GetKind_SU2() == SU2_DOT)) { Tot_Iter = System.Solve_b(StiffMatrix, LinSysRes, LinSysSol, geometry, config); - } su2double Residual = System.GetResidual(); diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 56141a8af063..aaf0e2eaf12a 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -90,7 +90,7 @@ template void CSysMatrix::Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar, unsigned short neqn, bool EdgeConnect, CGeometry *geometry, - CConfig *config, bool needTranspPtr) { + const CConfig *config, bool needTranspPtr) { assert(omp_get_thread_num()==0 && "Only the master thread is allowed to initialize the matrix."); @@ -107,15 +107,21 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi } /*--- Application of this matrix, FVM or FEM. ---*/ - auto type = EdgeConnect? ConnectivityType::FiniteVolume : ConnectivityType::FiniteElement; + const auto type = EdgeConnect? ConnectivityType::FiniteVolume : ConnectivityType::FiniteElement; - /*--- Types of preconditioner the matrix will be asked to build. ---*/ - unsigned short sol_prec = config->GetKind_Linear_Solver_Prec(); - unsigned short def_prec = config->GetKind_Deform_Linear_Solver_Prec(); - unsigned short adj_prec = config->GetKind_DiscAdj_Linear_Prec(); - bool adjoint = config->GetDiscrete_Adjoint(); + /*--- Type of preconditioner the matrix will be asked to build. ---*/ + auto prec = config->GetKind_Linear_Solver_Prec(); - bool ilu_needed = (sol_prec==ILU) || (def_prec==ILU) || (adjoint && (adj_prec==ILU)); + if (!EdgeConnect && !config->GetStructuralProblem()) { + /*--- FEM-type connectivity in non-structural context implies mesh deformation. ---*/ + prec = config->GetKind_Deform_Linear_Solver_Prec(); + } + else if (config->GetDiscrete_Adjoint() && (prec!=ILU)) { + /*--- Else "upgrade" primal solver settings. ---*/ + prec = config->GetKind_DiscAdj_Linear_Prec(); + } + const bool ilu_needed = (prec==ILU); + const bool diag_needed = ilu_needed || (prec==JACOBI) || (prec==LINELET); /*--- Basic dimensions. ---*/ nVar = nvar; @@ -166,9 +172,7 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi ALLOC_AND_INIT(ILU_matrix, nnz_ilu*nVar*nEqn) } - if (ilu_needed || (sol_prec==JACOBI) || (sol_prec==LINELET) || - (adjoint && (adj_prec==JACOBI)) || (def_prec==JACOBI)) - { + if (diag_needed) { ALLOC_AND_INIT(invM, nPointDomain*nVar*nEqn); } #undef ALLOC_AND_INIT @@ -1277,6 +1281,7 @@ void CSysMatrix::ComputeLineletPreconditioner(const CSysVector void CSysMatrix::ComputeResidual(const CSysVector & sol, const CSysVector & f, CSysVector & res) const { + SU2_OMP_BARRIER SU2_OMP_FOR_DYN(omp_heavy_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { ScalarType aux_vec[MAXNVAR]; diff --git a/Common/src/meson.build b/Common/src/meson.build index e1d65e7702c1..e311708e53d7 100644 --- a/Common/src/meson.build +++ b/Common/src/meson.build @@ -2,7 +2,7 @@ common_include = include_directories('./') common_src =files(['graph_coloring_structure.cpp', 'CConfig.cpp', 'blas_structure.cpp', - 'ad_structure.cpp', + 'basic_types/ad_structure.cpp', 'grid_movement_structure.cpp', 'grid_adaptation_structure.cpp', 'wall_model.cpp', diff --git a/SU2_CFD/include/solvers/CNSSolver.hpp b/SU2_CFD/include/solvers/CNSSolver.hpp index 060595e1f94d..0771a23c6940 100644 --- a/SU2_CFD/include/solvers/CNSSolver.hpp +++ b/SU2_CFD/include/solvers/CNSSolver.hpp @@ -78,10 +78,44 @@ class CNSSolver final : public CEulerSolver { unsigned long SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output) override; -protected: + /*! + * \brief Common code for wall boundaries, add the residual and Jacobian + * contributions due to grid motion associated with a particular boundary point. + */ + void AddDynamicGridResidualContribution(unsigned long iPoint, + unsigned long Point_Normal, + CGeometry* geometry, + const su2double* UnitNormal, + su2double Area, + const su2double* GridVel, + su2double** Jacobian_i, + su2double& Res_Conv, + su2double& Res_Visc) const; -public: + /*! + * \brief Get the wall temperature at a given vertex of a given marker for CHT problems. + */ + su2double GetCHTWallTemperature(const CConfig* config, + unsigned short val_marker, + unsigned long iVertex, + su2double thermal_conductivity, + su2double dist_ij, + su2double There, + su2double Temperature_Ref) const; + /*! + * \brief Generic implementation of the isothermal wall also covering CHT cases, + * for which the wall temperature is given by GetCHTWallTemperature. + */ + void BC_Isothermal_Wall_Generic(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker, + bool cht_mode = false); + +public: /*! * \brief Constructor of the class. */ diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index c627220b9208..9a8cc80fc32d 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -118,8 +118,6 @@ class CSolver { **Jacobian_jj; /*!< \brief Auxiliary matrices for storing point to point Jacobians. */ su2double *iPoint_UndLapl, /*!< \brief Auxiliary variable for the undivided Laplacians. */ *jPoint_UndLapl; /*!< \brief Auxiliary variable for the undivided Laplacians. */ - su2double **Smatrix, /*!< \brief Auxiliary structure for computing gradients by least-squares */ - **Cvector; /*!< \brief Auxiliary structure for computing gradients by least-squares */ int *Restart_Vars; /*!< \brief Auxiliary structure for holding the number of variables and points in a restart. */ int Restart_ExtIter; /*!< \brief Auxiliary structure for holding the external iteration offset from a restart. */ diff --git a/SU2_CFD/include/transport_model.hpp b/SU2_CFD/include/transport_model.hpp index 1c1a77ee039a..b27450f37de5 100644 --- a/SU2_CFD/include/transport_model.hpp +++ b/SU2_CFD/include/transport_model.hpp @@ -6,7 +6,7 @@ * * SU2 Project Website: https://su2code.github.io * - * The SU2 Project is maintained by the SU2 Foundation + * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md) @@ -41,7 +41,7 @@ #include "stdio.h" #include "math.h" -#include "../../Common/include/datatype_structure.hpp" +#include "../../Common/include/basic_types/datatype_structure.hpp" using namespace std; @@ -105,24 +105,24 @@ class CViscosityModel { * \version 1.0 */ class CConstantViscosity : public CViscosityModel { - + public: - + /*! * \brief Constructor of the class. */ CConstantViscosity(void); - + /*! * \brief Constructor of the class. */ CConstantViscosity(su2double mu_const); - + /*! * \brief Destructor of the class. */ ~CConstantViscosity(void) override; - + }; /*! @@ -137,34 +137,34 @@ class CSutherland : public CViscosityModel { su2double Mu_ref, /*!< \brief Internal Energy. */ T_ref, /*!< \brief DpDd_e. */ S; /*!< \brief DpDe_d. */ - + public: - + /*! * \brief Constructor of the class. */ CSutherland(void); - + /*! * \brief Constructor of the class. */ CSutherland(su2double mu_ref, su2double t_ref, su2double s); - + /*! * \brief Destructor of the class. */ ~CSutherland(void) override; - + /*! * \brief Set Viscosity. */ void SetViscosity(su2double T, su2double rho) override; - + /*! * \brief Set Viscosity Derivatives. */ void SetDerViscosity(su2double T, su2double rho) override; - + }; /*! @@ -178,27 +178,27 @@ class CPolynomialViscosity : public CViscosityModel { su2double *b; /*!< \brief Polynomial coefficients for viscosity as a function of temperature. */ public: - + /*! * \brief Constructor of the class. */ CPolynomialViscosity(void); - + /*! * \brief Constructor of the class. */ CPolynomialViscosity(unsigned short val_nCoeffs, const su2double* val_b); - + /*! * \brief Destructor of the class. */ ~CPolynomialViscosity(void) override; - + /*! * \brief Set Viscosity. */ void SetViscosity(su2double T, su2double rho) override; - + }; /*! @@ -286,33 +286,33 @@ class CConstantConductivity : public CConductivityModel { * \version 1.0 */ class CConstantConductivityRANS : public CConductivityModel { - + protected: su2double Kt_Lam; /*!< \brief Constant laminar conductivity. */ su2double Prandtl_Turb; /*!< \brief Turbulent Prandtl number. */ - + public: - + /*! * \brief Constructor of the class. */ CConstantConductivityRANS(void); - + /*! * \brief Constructor of the class. */ CConstantConductivityRANS(su2double kt_const, su2double pr_turb); - + /*! * \brief Destructor of the class. */ ~CConstantConductivityRANS(void) override; - + /*! * \brief Set effective thermal conductivity. */ void SetConductivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp) override; - + }; /*! @@ -401,29 +401,29 @@ class CPolynomialConductivity : public CConductivityModel { protected: unsigned short nPolyCoeffs; /*!< \brief Number of coefficients in the temperature polynomial. */ su2double *b; /*!< \brief Polynomial coefficients for thermal conductivity as a function of temperature. */ - + public: - + /*! * \brief Constructor of the class. */ CPolynomialConductivity(void); - + /*! * \brief Destructor of the class. */ ~CPolynomialConductivity(void) override; - + /*! * \brief Constructor of the class. */ CPolynomialConductivity(unsigned short val_nCoeffs, const su2double* val_b); - + /*! * \brief Set Thermal conductivity. */ void SetConductivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp) override; - + }; /*! @@ -438,27 +438,27 @@ class CPolynomialConductivityRANS : public CConductivityModel { su2double Prandtl_Turb; /*!< \brief Turbulent Prandtl number. */ public: - + /*! * \brief Constructor of the class. */ CPolynomialConductivityRANS(void); - + /*! * \brief Destructor of the class. */ ~CPolynomialConductivityRANS(void) override; - + /*! * \brief Constructor of the class. */ CPolynomialConductivityRANS(unsigned short val_nCoeffs, const su2double* val_b, su2double pr_turb); - + /*! * \brief Set Thermal conductivity. */ void SetConductivity(su2double T, su2double rho, su2double mu_lam, su2double mu_turb, su2double cp) override; - + }; #include "transport_model.inl" diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index b51b48313588..acc7f9e597a6 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -76,21 +76,21 @@ void CIntegration::Space_Integration(CGeometry *geometry, if (dual_time) solver_container[MainSolver]->SetResidual_DualTime(geometry, solver_container, config, iRKStep, iMesh, RunTime_EqSystem); - /// TODO: No boundary condition supports hybrid parallelism yet, master thread does all the work. + /*--- Pick convective and viscous numerics objects for the current thread. ---*/ - SU2_OMP_MASTER - { + CNumerics* conv_bound_numerics = numerics[CONV_BOUND_TERM + omp_get_thread_num()*MAX_TERMS]; + CNumerics* visc_bound_numerics = numerics[VISC_BOUND_TERM + omp_get_thread_num()*MAX_TERMS]; /*--- Boundary conditions that depend on other boundaries (they require MPI sincronization)---*/ - solver_container[MainSolver]->BC_Fluid_Interface(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config); + solver_container[MainSolver]->BC_Fluid_Interface(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config); /*--- Compute Fourier Transformations for markers where NRBC_BOUNDARY is applied---*/ if (config->GetBoolGiles() && config->GetSpatialFourier()){ - solver_container[MainSolver]->PreprocessBC_Giles(geometry, config, numerics[CONV_BOUND_TERM], INFLOW); + solver_container[MainSolver]->PreprocessBC_Giles(geometry, config, conv_bound_numerics, INFLOW); - solver_container[MainSolver]->PreprocessBC_Giles(geometry, config, numerics[CONV_BOUND_TERM], OUTFLOW); + solver_container[MainSolver]->PreprocessBC_Giles(geometry, config, conv_bound_numerics, OUTFLOW); } /*--- Weak boundary conditions ---*/ @@ -99,54 +99,54 @@ void CIntegration::Space_Integration(CGeometry *geometry, KindBC = config->GetMarker_All_KindBC(iMarker); switch (KindBC) { case EULER_WALL: - solver_container[MainSolver]->BC_Euler_Wall(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Euler_Wall(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case ACTDISK_INLET: - solver_container[MainSolver]->BC_ActDisk_Inlet(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_ActDisk_Inlet(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case ENGINE_INFLOW: - solver_container[MainSolver]->BC_Engine_Inflow(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Engine_Inflow(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case INLET_FLOW: - solver_container[MainSolver]->BC_Inlet(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Inlet(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case ACTDISK_OUTLET: - solver_container[MainSolver]->BC_ActDisk_Outlet(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_ActDisk_Outlet(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case ENGINE_EXHAUST: - solver_container[MainSolver]->BC_Engine_Exhaust(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Engine_Exhaust(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case SUPERSONIC_INLET: - solver_container[MainSolver]->BC_Supersonic_Inlet(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Supersonic_Inlet(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case OUTLET_FLOW: - solver_container[MainSolver]->BC_Outlet(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Outlet(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case SUPERSONIC_OUTLET: - solver_container[MainSolver]->BC_Supersonic_Outlet(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Supersonic_Outlet(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case GILES_BOUNDARY: - solver_container[MainSolver]->BC_Giles(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Giles(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case RIEMANN_BOUNDARY: if (config->GetBoolTurbomachinery()){ - solver_container[MainSolver]->BC_TurboRiemann(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_TurboRiemann(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); } else{ - solver_container[MainSolver]->BC_Riemann(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Riemann(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); } break; case FAR_FIELD: - solver_container[MainSolver]->BC_Far_Field(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Far_Field(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case SYMMETRY_PLANE: - solver_container[MainSolver]->BC_Sym_Plane(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Sym_Plane(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case ELECTRODE_BOUNDARY: - solver_container[MainSolver]->BC_Electrode(geometry, solver_container, numerics[CONV_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Electrode(geometry, solver_container, conv_bound_numerics, config, iMarker); break; case DIELEC_BOUNDARY: - solver_container[MainSolver]->BC_Dielec(geometry, solver_container, numerics[CONV_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Dielec(geometry, solver_container, conv_bound_numerics, config, iMarker); break; } } @@ -156,20 +156,20 @@ void CIntegration::Space_Integration(CGeometry *geometry, for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) switch (config->GetMarker_All_KindBC(iMarker)) { case ISOTHERMAL: - solver_container[MainSolver]->BC_Isothermal_Wall(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Isothermal_Wall(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case HEAT_FLUX: - solver_container[MainSolver]->BC_HeatFlux_Wall(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_HeatFlux_Wall(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case CUSTOM_BOUNDARY: - solver_container[MainSolver]->BC_Custom(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_Custom(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); break; case CHT_WALL_INTERFACE: if ((MainSolver == HEAT_SOL) || ((MainSolver == FLOW_SOL) && ((config->GetKind_Regime() == COMPRESSIBLE) || config->GetEnergy_Equation()))) { - solver_container[MainSolver]->BC_ConjugateHeat_Interface(geometry, solver_container, numerics[CONV_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_ConjugateHeat_Interface(geometry, solver_container, conv_bound_numerics, config, iMarker); } else { - solver_container[MainSolver]->BC_HeatFlux_Wall(geometry, solver_container, numerics[CONV_BOUND_TERM], numerics[VISC_BOUND_TERM], config, iMarker); + solver_container[MainSolver]->BC_HeatFlux_Wall(geometry, solver_container, conv_bound_numerics, visc_bound_numerics, config, iMarker); } break; } @@ -180,12 +180,9 @@ void CIntegration::Space_Integration(CGeometry *geometry, accumulated corectly during the communications. ---*/ if (config->GetnMarker_Periodic() > 0) { - solver_container[MainSolver]->BC_Periodic(geometry, solver_container, numerics[CONV_BOUND_TERM], config); + solver_container[MainSolver]->BC_Periodic(geometry, solver_container, conv_bound_numerics, config); } - } // end SU2_OMP_MASTER - SU2_OMP_BARRIER - } void CIntegration::Time_Integration(CGeometry *geometry, CSolver **solver_container, CConfig *config, diff --git a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp index 6b60823d06e9..465358a22530 100644 --- a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp @@ -191,18 +191,6 @@ CAdjEulerSolver::CAdjEulerSolver(CGeometry *geometry, CConfig *config, unsigned cout << "Explicit scheme. No Jacobian structure (Adjoint Euler). MG level: " << iMesh <<"." << endl; } - /*--- Computation of gradients by least squares ---*/ - if (config->GetLeastSquaresRequired()) { - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - /*--- c vector := transpose(WA)*(Wb) ---*/ - Cvector = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - /*--- Sensitivity definition and coefficient in all the markers ---*/ CSensitivity = new su2double* [nMarker]; for (iMarker = 0; iMarker < nMarker; iMarker++) { diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp index 742c51230b02..332bed424799 100644 --- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp @@ -151,18 +151,6 @@ CAdjNSSolver::CAdjNSSolver(CGeometry *geometry, CConfig *config, unsigned short cout << "Explicit scheme. No Jacobian structure (Adjoint N-S). MG level: " << iMesh <<"." << endl; } - /*--- Array structures for computation of gradients by least squares ---*/ - if (config->GetLeastSquaresRequired()) { - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - /*--- c vector := transpose(WA)*(Wb) ---*/ - Cvector = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - /*--- Sensitivity definition and coefficient on all markers ---*/ CSensitivity = new su2double* [nMarker]; for (iMarker=0; iMarkerGetLeastSquaresRequired()) { - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - /*--- c vector := transpose(WA)*(Wb) ---*/ - Cvector = new su2double* [nVar+1]; - for (iVar = 0; iVar < nVar+1; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - /*--- Far-Field values and initizalization ---*/ bool restart = config->GetRestart(); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 7a9706cee858..c22ec7d41f8a 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -169,14 +169,8 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, /*--- Define some auxiliar vector related with the residual ---*/ - Residual = new su2double[nVar](); Residual_RMS = new su2double[nVar](); Residual_Max = new su2double[nVar](); - Residual_i = new su2double[nVar](); - Residual_j = new su2double[nVar](); - Res_Conv = new su2double[nVar](); - Res_Visc = new su2double[nVar](); - Res_Sour = new su2double[nVar](); /*--- Define some structures for locating max residuals ---*/ @@ -186,14 +180,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Point_Max_Coord[iVar] = new su2double[nDim](); } - /*--- Define some auxiliary vectors related to the solution ---*/ - - Solution = new su2double[nVar](); - - /*--- Define some auxiliary vectors related to the geometry ---*/ - - Vector = new su2double[nDim](); - /*--- Define some auxiliar vector related with the undivided lapalacian computation ---*/ if (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) { @@ -263,13 +249,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, if (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT) { - Jacobian_i = new su2double* [nVar]; - Jacobian_j = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double [nVar]; - Jacobian_j[iVar] = new su2double [nVar]; - } - if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (" << description << "). MG level: " << iMesh <<"." << endl; @@ -288,20 +267,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, cout << "Explicit scheme. No Jacobian structure (" << description << "). MG level: " << iMesh <<"." << endl; } - /*--- Define some auxiliary vectors for computing flow variable - gradients by least squares, S matrix := inv(R)*traspose(inv(R)), - c vector := transpose(WA)*(Wb) ---*/ - - if (config->GetLeastSquaresRequired()) { - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - - Cvector = new su2double* [nPrimVarGrad]; - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - /*--- Allocates a 2D array with variable "outer" sizes and init to 0. ---*/ auto Alloc2D = [](unsigned long M, const unsigned long* N, su2double**& X) { @@ -491,6 +456,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, /*--- Use the values at the infinity ---*/ + su2double Solution[MAXNVAR] = {0.0}; if ((Pressure < 0.0) || (Density < 0.0) || (Temperature < 0.0)) { Solution[0] = Density_Inf; for (iDim = 0; iDim < nDim; iDim++) @@ -3386,9 +3352,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Include the residual contribution from GCL due to the static mesh movement that is set for rotating frame. ---*/ - SU2_OMP_MASTER SetRotatingFrame_GCL(geometry, config); - SU2_OMP_BARRIER /*--- Loop over all points ---*/ SU2_OMP_FOR_DYN(omp_chunk_size) @@ -4901,6 +4865,7 @@ void CEulerSolver::GetPower_Properties(CGeometry *geometry, CConfig *config, uns su2double DeltaPress = 0.0, DeltaTemp = 0.0, TotalPressRatio = 0.0, TotalTempRatio = 0.0, StaticPressRatio = 0.0, StaticTempRatio = 0.0, NetThrust = 0.0, GrossThrust = 0.0, Power = 0.0, MassFlow = 0.0, Mach = 0.0, Force = 0.0; bool ReverseFlow, Engine = false, Pair = true; + su2double Vector[MAXNDIM] = {0.0}; su2double Gas_Constant = config->GetGas_ConstantND(); su2double Cp = Gas_Constant*Gamma / (Gamma-1.0); @@ -5795,6 +5760,7 @@ void CEulerSolver::SetActDisk_BCThrust(CGeometry *geometry, CSolver **solver_con su2double Target_Force, Force, Target_Power, Power, NetThrust, BCThrust_old, Initial_BCThrust; bool ActDisk_Info; su2double MyBCThrust, BCThrust_Init; + su2double Vector[MAXNDIM] = {0.0}; su2double dNetThrust_dBCThrust = config->GetdNetThrust_dBCThrust(); unsigned short Kind_ActDisk = config->GetKind_ActDisk(); @@ -6747,21 +6713,18 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, unsigned short iDim, iVar; unsigned long iVertex, iPoint; - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT), - viscous = config->GetViscous(); + bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + bool viscous = config->GetViscous(); + bool preprocessed = false; /*--- Allocation of variables necessary for convective fluxes. ---*/ - su2double Area, ProjVelocity_i, - *V_reflected, - *V_domain, - *Normal = new su2double[nDim], - *UnitNormal = new su2double[nDim]; + su2double Area, ProjVelocity_i, *V_reflected, *V_domain, + Normal[MAXNDIM] = {0.0}, UnitNormal[MAXNDIM] = {0.0}; /*--- Allocation of variables necessary for viscous fluxes. ---*/ su2double ProjGradient, ProjNormVelGrad, ProjTangVelGrad, TangentialNorm, - *Tangential = new su2double[nDim], - *GradNormVel = new su2double[nDim], - *GradTangVel = new su2double[nDim]; + Tangential[MAXNDIM] = {0.0}, GradNormVel[MAXNDIM] = {0.0}, + GradTangVel[MAXNDIM] = {0.0}; /*--- Allocation of primitive gradient arrays for viscous fluxes. ---*/ su2double **Grad_Reflected = new su2double*[nPrimVarGrad]; @@ -6769,9 +6732,11 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, Grad_Reflected[iVar] = new su2double[nDim]; /*--- Loop over all the vertices on this boundary marker. ---*/ + + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - if (iVertex == 0 || + if (!preprocessed || geometry->bound_is_straight[val_marker] != true) { /*----------------------------------------------------------------------------------------------*/ @@ -6787,6 +6752,8 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- penalty). ---*/ /*----------------------------------------------------------------------------------------------*/ + preprocessed = true; + /*--- Normal vector for a random vertex (zero) on this marker (negate for outward convention). ---*/ geometry->vertex[val_marker][iVertex]->GetNormal(Normal); for (iDim = 0; iDim < nDim; iDim++) @@ -6982,7 +6949,7 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), Grad_Reflected); /*--- Turbulent kinetic energy. ---*/ - if (config->GetKind_Turb_Model() == SST) + if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); @@ -7000,12 +6967,6 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, }//for iVertex /*--- Free locally allocated memory ---*/ - delete [] Normal; - delete [] UnitNormal; - delete [] Tangential; - delete [] GradNormVel; - delete [] GradTangVel; - for (iVar = 0; iVar < nPrimVarGrad; iVar++) delete [] Grad_Reflected[iVar]; delete [] Grad_Reflected; @@ -7039,6 +7000,7 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, /*--- Loop over all the vertices on this boundary marker ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); @@ -7282,8 +7244,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, su2double *Velocity_e, Velocity2_e, VelMag_e, Enthalpy_e, Entropy_e, Energy_e = 0.0, StaticEnthalpy_e, StaticEnergy_e, Density_e = 0.0, Pressure_e; su2double *Velocity_i, Velocity2_i, Enthalpy_i, Energy_i, StaticEnergy_i, Density_i, Kappa_i, Chi_i, Pressure_i, SoundSpeed_i; su2double ProjVelocity_i; - su2double **P_Tensor, **invP_Tensor, *Lambda_i, **Jacobian_b, **DubDu, *dw, *u_e, *u_i, *u_b; - su2double *gridVel; + su2double **P_Tensor, **invP_Tensor, *Lambda_i, **Jacobian_b, **Jacobian_i, **DubDu, *dw, *u_e, *u_i, *u_b; + su2double *gridVel, *Residual; su2double *V_boundary, *V_domain, *S_boundary, *S_domain; bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); @@ -7305,17 +7267,23 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, u_b = new su2double[nVar]; dw = new su2double[nVar]; + Residual = new su2double[nVar]; + S_boundary = new su2double[8]; P_Tensor = new su2double*[nVar]; invP_Tensor = new su2double*[nVar]; + Jacobian_i = new su2double*[nVar]; for (iVar = 0; iVar < nVar; iVar++) { P_Tensor[iVar] = new su2double[nVar]; invP_Tensor[iVar] = new su2double[nVar]; + Jacobian_i[iVar] = new su2double[nVar]; } /*--- Loop over all the vertices on this boundary marker ---*/ + + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { V_boundary= GetCharacPrimVar(val_marker, iVertex); @@ -7763,14 +7731,17 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, delete [] u_b; delete [] dw; + delete [] Residual; for (iVar = 0; iVar < nVar; iVar++) { delete [] P_Tensor[iVar]; delete [] invP_Tensor[iVar]; + delete [] Jacobian_i[iVar]; } delete [] P_Tensor; delete [] invP_Tensor; + delete [] Jacobian_i; } @@ -7786,8 +7757,8 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain su2double *Velocity_e, Velocity2_e, Enthalpy_e, Entropy_e, Energy_e = 0.0, StaticEnthalpy_e, StaticEnergy_e, Density_e = 0.0, Pressure_e; su2double *Velocity_i, Velocity2_i, Enthalpy_i, Energy_i, StaticEnergy_i, Density_i, Kappa_i, Chi_i, Pressure_i, SoundSpeed_i; su2double ProjVelocity_i; - su2double **P_Tensor, **invP_Tensor, *Lambda_i, **Jacobian_b, **DubDu, *dw, *u_e, *u_i, *u_b; - su2double *gridVel; + su2double **P_Tensor, **invP_Tensor, *Lambda_i, **Jacobian_b, **Jacobian_i, **DubDu, *dw, *u_e, *u_i, *u_b; + su2double *gridVel, *Residual; su2double *V_boundary, *V_domain, *S_boundary, *S_domain; su2double AverageEnthalpy, AverageEntropy; unsigned short iZone = config->GetiZone(); @@ -7814,18 +7785,24 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain u_b = new su2double[nVar]; dw = new su2double[nVar]; + Residual = new su2double[nVar]; + S_boundary = new su2double[8]; P_Tensor = new su2double*[nVar]; invP_Tensor = new su2double*[nVar]; + Jacobian_i = new su2double*[nVar]; for (iVar = 0; iVar < nVar; iVar++) { P_Tensor[iVar] = new su2double[nVar]; invP_Tensor[iVar] = new su2double[nVar]; + Jacobian_i[iVar] = new su2double[nVar]; } /*--- Loop over all the vertices on this boundary marker ---*/ for (iSpan= 0; iSpan < nSpanWiseSections; iSpan++){ + + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { /*--- using the other vertex information for retrieving some information ---*/ @@ -8270,14 +8247,17 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain delete [] u_b; delete [] dw; + delete [] Residual; for (iVar = 0; iVar < nVar; iVar++) { delete [] P_Tensor[iVar]; delete [] invP_Tensor[iVar]; + delete [] Jacobian_i[iVar]; } delete [] P_Tensor; delete [] invP_Tensor; + delete [] Jacobian_i; } @@ -8749,6 +8729,7 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu /*--- Loop over all the vertices on this boundary marker ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { /*--- using the other vertex information for retrieving some information ---*/ @@ -9198,6 +9179,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, /*--- Loop over all the vertices on this boundary marker ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Allocate the value at the inlet ---*/ @@ -9493,6 +9475,8 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, su2double *Normal = new su2double[nDim]; /*--- Loop over all the vertices on this boundary marker ---*/ + + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Allocate the value at the outlet ---*/ @@ -9693,6 +9677,7 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con /*--- Loop over all the vertices on this boundary marker ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Allocate the value at the outlet ---*/ @@ -9810,6 +9795,7 @@ void CEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solver_co /*--- Loop over all the vertices on this boundary marker ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); @@ -9992,6 +9978,7 @@ void CEulerSolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_contai /*--- Loop over all the vertices on this boundary marker ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Allocate the value at the outlet ---*/ @@ -10180,6 +10167,7 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta /*--- Loop over all the vertices on this boundary marker ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Allocate the value at the exhaust ---*/ @@ -10403,6 +10391,10 @@ void CEulerSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_cont su2double PrimVar_i[MAXNVAR] = {0.0}; su2double PrimVar_j[MAXNVAR] = {0.0}; su2double Secondary_j[MAXNVAR] = {0.0}; + su2double Residual[MAXNVAR] = {0.0}; + su2double **Jacobian_i = new su2double* [nVar]; + for (iVar = 0; iVar < nVar; iVar++) + Jacobian_i[iVar] = new su2double [nVar]; su2double weight; su2double P_static, rho_static; @@ -10411,6 +10403,7 @@ void CEulerSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_cont if (config->GetMarker_All_KindBC(iMarker) == FLUID_INTERFACE) { + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); @@ -10555,6 +10548,10 @@ void CEulerSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_cont } } + for (iVar = 0; iVar < nVar; iVar++) + delete [] Jacobian_i[iVar]; + delete [] Jacobian_i; + } void CEulerSolver::BC_Interface_Boundary(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, @@ -10572,6 +10569,7 @@ void CEulerSolver::BC_Interface_Boundary(CGeometry *geometry, CSolver **solver_c /*--- Do the send process, by the moment we are sending each node individually, this must be changed ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); @@ -10634,6 +10632,7 @@ void CEulerSolver::BC_NearField_Boundary(CGeometry *geometry, CSolver **solver_c /*--- Do the send process, by the moment we are sending each node individually, this must be changed ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); @@ -10725,6 +10724,7 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C /*--- Loop over all the vertices on this boundary marker ---*/ + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); @@ -11117,10 +11117,12 @@ void CEulerSolver::BC_Periodic(CGeometry *geometry, CSolver **solver_container, accumulated correctly during the communications. For implicit calculations, the Jacobians and linear system are also correctly adjusted here. ---*/ + SU2_OMP_MASTER for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_RESIDUAL); CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_RESIDUAL); } + SU2_OMP_BARRIER } @@ -11147,6 +11149,7 @@ void CEulerSolver::BC_Custom(CGeometry *geometry, /*--- Loop over all the vertices on this boundary marker ---*/ + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { /*--- Get the point index for the current node. ---*/ @@ -11163,6 +11166,7 @@ void CEulerSolver::BC_Custom(CGeometry *geometry, /*--- Get the conservative state from the verification solution. ---*/ + su2double Solution[MAXNVAR] = {0.0}; VerificationSolution->GetBCState(coor, time, Solution); /*--- For verification cases, we will apply a strong Dirichlet diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 75eefd05ab95..aab05b51cab9 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1952,13 +1952,16 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CSysVector LinSysAux(nPoint, nPointDomain, nVar, nullptr); #endif +#if defined(CODI_REVERSE_TYPE) || defined(USE_MIXED_PRECISION) + /*--- We need temporaries to interface with the passive matrix. ---*/ + CSysVector sol, res; +#endif + SU2_OMP_PARALLEL { #if !defined(CODI_REVERSE_TYPE) && !defined(USE_MIXED_PRECISION) Jacobian.ComputeResidual(LinSysSol, LinSysRes, LinSysAux); #else - /*--- We need temporaries to interface with the passive matrix. ---*/ - CSysVector sol, res; sol.PassiveCopy(LinSysSol); res.PassiveCopy(LinSysRes); Jacobian.ComputeResidual(sol, res, LinSysAux); diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index 418a4fdb2e26..53dcba53a07f 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -134,19 +134,6 @@ CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iM OutputHeadingNames = new string[nOutputVariables]; } - /*--- Computation of gradients by least squares ---*/ - - if (config->GetLeastSquaresRequired()) { - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - /*--- c vector := transpose(WA)*(Wb) ---*/ - Cvector = new su2double* [nVar+1]; - for (iVar = 0; iVar < nVar+1; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - HeatFlux_per_Marker = new su2double[nMarker]; AverageT_per_Marker = new su2double[nMarker]; Surface_Areas = new su2double[config->GetnMarker_HeatFlux()]; diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 53e523240386..4d0becdb1a77 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -74,7 +74,6 @@ CIncEulerSolver::CIncEulerSolver(void) : CSolver() { jPoint_UndLapl = nullptr; Primitive = nullptr; Primitive_i = nullptr; Primitive_j = nullptr; CharacPrimVar = nullptr; - Smatrix = nullptr; Cvector = nullptr; Preconditioner = nullptr; FluidModel = nullptr; @@ -183,7 +182,6 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned jPoint_UndLapl = nullptr; Primitive = nullptr; Primitive_i = nullptr; Primitive_j = nullptr; CharacPrimVar = nullptr; - Smatrix = nullptr; Cvector = nullptr; Preconditioner = nullptr; /*--- Fluid model pointer initialization ---*/ @@ -310,22 +308,6 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned if (rank == MASTER_NODE) cout << "Explicit scheme. No Jacobian structure (Euler). MG level: " << iMesh <<"." << endl; } - /*--- Define some auxiliary vectors for computing flow variable - gradients by least squares, S matrix := inv(R)*traspose(inv(R)), - c vector := transpose(WA)*(Wb) ---*/ - - if (config->GetLeastSquaresRequired()) { - - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - - Cvector = new su2double* [nPrimVarGrad]; - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - Cvector[iVar] = new su2double [nDim]; - - } - /*--- Store the value of the characteristic primitive variables at the boundaries ---*/ CharacPrimVar = new su2double** [nMarker]; diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 258c433cd512..e5aa616c3b40 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -252,20 +252,6 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short cout << "Explicit scheme. No Jacobian structure (Navier-Stokes). MG level: " << iMesh <<"." << endl; } - /*--- Define some auxiliary vectors for computing flow variable - gradients by least squares, S matrix := inv(R)*traspose(inv(R)), - c vector := transpose(WA)*(Wb) ---*/ - - if (config->GetLeastSquaresRequired()) { - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - - Cvector = new su2double* [nPrimVarGrad]; - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - /*--- Store the value of the characteristic primitive variables at the boundaries ---*/ CharacPrimVar = new su2double** [nMarker]; diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index a860bc4dc03b..e596a4574729 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -25,10 +25,11 @@ * License along with SU2. If not, see . */ - #include "../../include/solvers/CNSSolver.hpp" #include "../../include/variables/CNSVariable.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" + CNSSolver::CNSSolver(void) : CEulerSolver() { } @@ -256,9 +257,7 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C /*--- Compute the TauWall from the wall functions ---*/ if (wall_functions) { - SU2_OMP_MASTER SetTauWall_WF(geometry, solver_container, config); - SU2_OMP_BARRIER } } @@ -941,1120 +940,648 @@ void CNSSolver::Evaluate_ObjFunc(CConfig *config) { } +void CNSSolver::SetRoe_Dissipation(CGeometry *geometry, CConfig *config){ -void CNSSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - - unsigned short iDim, jDim, iVar, jVar; - unsigned long iVertex, iPoint, Point_Normal, total_index; - - su2double Wall_HeatFlux, dist_ij, *Coord_i, *Coord_j, theta2; - su2double thetax, thetay, thetaz, etax, etay, etaz, pix, piy, piz, factor; - su2double ProjGridVel, *GridVel, GridVel2, *Normal, Area, Pressure = 0.0; - su2double total_viscosity, div_vel, Density, tau_vel[3] = {0.0, 0.0, 0.0}, UnitNormal[3] = {0.0, 0.0, 0.0}; - su2double laminar_viscosity = 0.0, eddy_viscosity = 0.0, Grad_Vel[3][3] = {{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}}, - tau[3][3] = {{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}}; - su2double delta[3][3] = {{1.0, 0.0, 0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}}; - - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - - /*--- Identify the boundary by string name ---*/ - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - - /*--- Get the specified wall heat flux from config as well as the - wall function treatment.---*/ - - Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag)/config->GetHeat_Flux_Ref(); - -// Wall_Function = config->GetWallFunction_Treatment(Marker_Tag); -// if (Wall_Function != NO_WALL_FUNCTION) { -// SU2_MPI::Error("Wall function treament not implemented yet", CURRENT_FUNCTION); -// } - - /*--- Loop over all of the vertices on this boundary marker ---*/ - - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- If it is a customizable patch, retrieve the specified wall heat flux. ---*/ - - if (config->GetMarker_All_PyCustom(val_marker)) Wall_HeatFlux = geometry->GetCustomBoundaryHeatFlux(val_marker, iVertex); - - /*--- Compute dual-grid area and boundary normal ---*/ - - Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - - Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Area += Normal[iDim]*Normal[iDim]; - Area = sqrt (Area); - - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = -Normal[iDim]/Area; - - /*--- Initialize the convective & viscous residuals to zero ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Res_Conv[iVar] = 0.0; - Res_Visc[iVar] = 0.0; - } - - /*--- Store the corrected velocity at the wall which will - be zero (v = 0), unless there are moving walls (v = u_wall)---*/ - - if (dynamic_grid) { - GridVel = geometry->nodes->GetGridVel(iPoint); - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = GridVel[iDim]; - } else { - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = 0.0; - } - - /*--- Impose the value of the velocity as a strong boundary - condition (Dirichlet). Fix the velocity and remove any - contribution to the residual at this node. ---*/ - - nodes->SetVelocity_Old(iPoint,Vector); - - for (iDim = 0; iDim < nDim; iDim++) - LinSysRes.SetBlock_Zero(iPoint, iDim+1); - nodes->SetVel_ResTruncError_Zero(iPoint); - - /*--- Apply a weak boundary condition for the energy equation. - Compute the residual due to the prescribed heat flux. ---*/ - - Res_Visc[nDim+1] = Wall_HeatFlux * Area; - - /*--- If the wall is moving, there are additional residual contributions - due to pressure (p v_wall.n) and shear stress (tau.v_wall.n). ---*/ - - if (dynamic_grid) { - - /*--- Get the grid velocity at the current boundary node ---*/ - - GridVel = geometry->nodes->GetGridVel(iPoint); - ProjGridVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjGridVel += GridVel[iDim]*UnitNormal[iDim]*Area; - - /*--- Retrieve other primitive quantities and viscosities ---*/ - - Density = nodes->GetDensity(iPoint); - Pressure = nodes->GetPressure(iPoint); - laminar_viscosity = nodes->GetLaminarViscosity(iPoint); - eddy_viscosity = nodes->GetEddyViscosity(iPoint); - total_viscosity = laminar_viscosity + eddy_viscosity; - - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint,iDim+1, jDim); - } - } - - /*--- Divergence of the velocity ---*/ + const unsigned short kind_roe_dissipation = config->GetKind_RoeLowDiss(); - div_vel = 0.0; for (iDim = 0 ; iDim < nDim; iDim++) div_vel += Grad_Vel[iDim][iDim]; + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - /*--- Compute the viscous stress tensor ---*/ + if (kind_roe_dissipation == FD || kind_roe_dissipation == FD_DUCROS){ - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0; jDim < nDim; jDim++) { - tau[iDim][jDim] = total_viscosity*( Grad_Vel[jDim][iDim]+Grad_Vel[iDim][jDim] ) - - TWO3*total_viscosity*div_vel*delta[iDim][jDim]; - } - } + su2double wall_distance = geometry->nodes->GetWall_Distance(iPoint); - /*--- Dot product of the stress tensor with the grid velocity ---*/ + nodes->SetRoe_Dissipation_FD(iPoint, wall_distance); - for (iDim = 0 ; iDim < nDim; iDim++) { - tau_vel[iDim] = 0.0; - for (jDim = 0 ; jDim < nDim; jDim++) - tau_vel[iDim] += tau[iDim][jDim]*GridVel[jDim]; - } + } else if (kind_roe_dissipation == NTS || kind_roe_dissipation == NTS_DUCROS) { - /*--- Compute the convective and viscous residuals (energy eqn.) ---*/ + const su2double delta = geometry->nodes->GetMaxLength(iPoint); + assert(delta > 0 && "Delta must be initialized and non-negative"); + nodes->SetRoe_Dissipation_NTS(iPoint, delta, config->GetConst_DES()); + } + } - Res_Conv[nDim+1] = Pressure*ProjGridVel; - for (iDim = 0 ; iDim < nDim; iDim++) - Res_Visc[nDim+1] += tau_vel[iDim]*UnitNormal[iDim]*Area; +} - /*--- Implicit Jacobian contributions due to moving walls ---*/ +void CNSSolver::AddDynamicGridResidualContribution(unsigned long iPoint, unsigned long Point_Normal, + CGeometry* geometry, const su2double* UnitNormal, + su2double Area, const su2double* GridVel, + su2double** Jacobian_i, su2double& Res_Conv, + su2double& Res_Visc) const { - if (implicit) { + su2double ProjGridVel = Area * GeometryToolbox::DotProduct(nDim, GridVel, UnitNormal); - /*--- Jacobian contribution related to the pressure term ---*/ + /*--- Retrieve other primitive quantities and viscosities ---*/ - GridVel2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - GridVel2 += GridVel[iDim]*GridVel[iDim]; - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - Jacobian_i[nDim+1][0] = 0.5*(Gamma-1.0)*GridVel2*ProjGridVel; - for (jDim = 0; jDim < nDim; jDim++) - Jacobian_i[nDim+1][jDim+1] = -(Gamma-1.0)*GridVel[jDim]*ProjGridVel; - Jacobian_i[nDim+1][nDim+1] = (Gamma-1.0)*ProjGridVel; + su2double Density = nodes->GetDensity(iPoint); + su2double Pressure = nodes->GetPressure(iPoint); + su2double laminar_viscosity = nodes->GetLaminarViscosity(iPoint); + su2double eddy_viscosity = nodes->GetEddyViscosity(iPoint); + su2double total_viscosity = laminar_viscosity + eddy_viscosity; - /*--- Add the block to the Global Jacobian structure ---*/ + const auto Grad_Vel = &nodes->GetGradient_Primitive(iPoint)[1]; - Jacobian.AddBlock2Diag(iPoint, Jacobian_i); + /*--- Divergence of the velocity ---*/ - /*--- Now the Jacobian contribution related to the shear stress ---*/ + su2double div_vel = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) + div_vel += Grad_Vel[iDim][iDim]; - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; + /*--- Compute the viscous stress tensor ---*/ - /*--- Compute closest normal neighbor ---*/ + su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}; + for (auto iDim = 0u; iDim < nDim; iDim++) { + for (auto jDim = 0u; jDim < nDim; jDim++) { + tau[iDim][jDim] = total_viscosity * (Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim]); + } + tau[iDim][iDim] -= TWO3*total_viscosity*div_vel; + } - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + /*--- Dot product of the stress tensor with the grid velocity ---*/ - /*--- Get coordinates of i & nearest normal and compute distance ---*/ + su2double tau_vel[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + tau_vel[iDim] = GeometryToolbox::DotProduct(nDim, tau[iDim], GridVel); - Coord_i = geometry->nodes->GetCoord(iPoint); - Coord_j = geometry->nodes->GetCoord(Point_Normal); + /*--- Compute the convective and viscous residuals (energy eqn.) ---*/ - dist_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) - dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - dist_ij = sqrt(dist_ij); + Res_Conv += Pressure*ProjGridVel; + Res_Visc += GeometryToolbox::DotProduct(nDim, tau_vel, UnitNormal) * Area; - theta2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - theta2 += UnitNormal[iDim]*UnitNormal[iDim]; + /*--- Implicit Jacobian contributions due to moving walls ---*/ - factor = total_viscosity*Area/(Density*dist_ij); + if (Jacobian_i != nullptr) { - if (nDim == 2) { - thetax = theta2 + UnitNormal[0]*UnitNormal[0]/3.0; - thetay = theta2 + UnitNormal[1]*UnitNormal[1]/3.0; + /*--- Jacobian contribution related to the pressure term ---*/ - etaz = UnitNormal[0]*UnitNormal[1]/3.0; + su2double GridVel2 = GeometryToolbox::SquaredNorm(nDim, GridVel); - pix = GridVel[0]*thetax + GridVel[1]*etaz; - piy = GridVel[0]*etaz + GridVel[1]*thetay; + Jacobian_i[nDim+1][0] += 0.5*(Gamma-1.0)*GridVel2*ProjGridVel; - Jacobian_i[nDim+1][0] -= factor*(-pix*GridVel[0]+piy*GridVel[1]); - Jacobian_i[nDim+1][1] -= factor*pix; - Jacobian_i[nDim+1][2] -= factor*piy; - } else { - thetax = theta2 + UnitNormal[0]*UnitNormal[0]/3.0; - thetay = theta2 + UnitNormal[1]*UnitNormal[1]/3.0; - thetaz = theta2 + UnitNormal[2]*UnitNormal[2]/3.0; + for (auto jDim = 0u; jDim < nDim; jDim++) + Jacobian_i[nDim+1][jDim+1] += -(Gamma-1.0)*GridVel[jDim]*ProjGridVel; - etaz = UnitNormal[0]*UnitNormal[1]/3.0; - etax = UnitNormal[1]*UnitNormal[2]/3.0; - etay = UnitNormal[0]*UnitNormal[2]/3.0; + Jacobian_i[nDim+1][nDim+1] += (Gamma-1.0)*ProjGridVel; - pix = GridVel[0]*thetax + GridVel[1]*etaz + GridVel[2]*etay; - piy = GridVel[0]*etaz + GridVel[1]*thetay + GridVel[2]*etax; - piz = GridVel[0]*etay + GridVel[1]*etax + GridVel[2]*thetaz; + /*--- Now the Jacobian contribution related to the shear stress ---*/ - Jacobian_i[nDim+1][0] -= factor*(-pix*GridVel[0]+piy*GridVel[1]+piz*GridVel[2]); - Jacobian_i[nDim+1][1] -= factor*pix; - Jacobian_i[nDim+1][2] -= factor*piy; - Jacobian_i[nDim+1][3] -= factor*piz; - } + /*--- Get coordinates of i & nearest normal and compute distance ---*/ - /*--- Subtract the block from the Global Jacobian structure ---*/ + const auto Coord_i = geometry->nodes->GetCoord(iPoint); + const auto Coord_j = geometry->nodes->GetCoord(Point_Normal); - Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); + su2double dist_ij = GeometryToolbox::Distance(nDim, Coord_i, Coord_j); - } - } + const su2double theta2 = 1.0; - /*--- Convective contribution to the residual at the wall ---*/ + su2double factor = total_viscosity*Area/(Density*dist_ij); - LinSysRes.AddBlock(iPoint, Res_Conv); + if (nDim == 2) { + su2double thetax = theta2 + UnitNormal[0]*UnitNormal[0]/3.0; + su2double thetay = theta2 + UnitNormal[1]*UnitNormal[1]/3.0; - /*--- Viscous contribution to the residual at the wall ---*/ + su2double etaz = UnitNormal[0]*UnitNormal[1]/3.0; - LinSysRes.SubtractBlock(iPoint, Res_Visc); - - /*--- Enforce the no-slip boundary condition in a strong way by - modifying the velocity-rows of the Jacobian (1 on the diagonal). ---*/ - - if (implicit) { - for (iVar = 1; iVar <= nDim; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } - } + su2double pix = GridVel[0]*thetax + GridVel[1]*etaz; + su2double piy = GridVel[0]*etaz + GridVel[1]*thetay; + Jacobian_i[nDim+1][0] += factor*(-pix*GridVel[0]+piy*GridVel[1]); + Jacobian_i[nDim+1][1] += factor*pix; + Jacobian_i[nDim+1][2] += factor*piy; + } + else { + su2double thetax = theta2 + UnitNormal[0]*UnitNormal[0]/3.0; + su2double thetay = theta2 + UnitNormal[1]*UnitNormal[1]/3.0; + su2double thetaz = theta2 + UnitNormal[2]*UnitNormal[2]/3.0; + + su2double etaz = UnitNormal[0]*UnitNormal[1]/3.0; + su2double etax = UnitNormal[1]*UnitNormal[2]/3.0; + su2double etay = UnitNormal[0]*UnitNormal[2]/3.0; + + su2double pix = GridVel[0]*thetax + GridVel[1]*etaz + GridVel[2]*etay; + su2double piy = GridVel[0]*etaz + GridVel[1]*thetay + GridVel[2]*etax; + su2double piz = GridVel[0]*etay + GridVel[1]*etax + GridVel[2]*thetaz; + + Jacobian_i[nDim+1][0] += factor*(-pix*GridVel[0]+piy*GridVel[1]+piz*GridVel[2]); + Jacobian_i[nDim+1][1] += factor*pix; + Jacobian_i[nDim+1][2] += factor*piy; + Jacobian_i[nDim+1][3] += factor*piz; } } } -void CNSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - - unsigned short iVar, jVar, iDim, jDim; - unsigned long iVertex, iPoint, Point_Normal, total_index; - - su2double *Normal, *Coord_i, *Coord_j, Area, dist_ij, theta2; - su2double Twall, dTdn, dTdrho, thermal_conductivity; - su2double thetax, thetay, thetaz, etax, etay, etaz, pix, piy, piz, factor; - su2double ProjGridVel, *GridVel, GridVel2, Pressure = 0.0, Density, Vel2; - su2double total_viscosity, div_vel, tau_vel[3] = {0.0,0.0,0.0}, UnitNormal[3] = {0.0,0.0,0.0}; - su2double laminar_viscosity, eddy_viscosity, Grad_Vel[3][3] = {{1.0, 0.0, 0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}}, - tau[3][3] = {{0.0, 0.0, 0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}}, delta[3][3] = {{1.0, 0.0, 0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}}; - - su2double Prandtl_Lam = config->GetPrandtl_Lam(); - su2double Prandtl_Turb = config->GetPrandtl_Turb(); - su2double Gas_Constant = config->GetGas_ConstantND(); - su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; - - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - - /*--- Identify the boundary ---*/ - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); +void CNSSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - /*--- Retrieve the specified wall temperature from config - as well as the wall function treatment.---*/ + /*--- Identify the boundary by string name and get the specified wall + heat flux from config as well as the wall function treatment. ---*/ - Twall = config->GetIsothermal_Temperature(Marker_Tag)/config->GetTemperature_Ref(); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); + su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag)/config->GetHeat_Flux_Ref(); // Wall_Function = config->GetWallFunction_Treatment(Marker_Tag); // if (Wall_Function != NO_WALL_FUNCTION) { // SU2_MPI::Error("Wall function treament not implemented yet", CURRENT_FUNCTION); // } - /*--- Loop over boundary points ---*/ - - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- If it is a customizable patch, retrieve the specified wall temperature. ---*/ - - if (config->GetMarker_All_PyCustom(val_marker)) Twall = geometry->GetCustomBoundaryTemperature(val_marker, iVertex); - - /*--- Compute dual-grid area and boundary normal ---*/ - - Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - - Area = 0.0; for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; Area = sqrt (Area); - - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = -Normal[iDim]/Area; - - /*--- Calculate useful quantities ---*/ - - theta2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - theta2 += UnitNormal[iDim]*UnitNormal[iDim]; - - /*--- Compute closest normal neighbor ---*/ - - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + /*--- Jacobian, initialized to zero if needed. ---*/ + su2double **Jacobian_i = nullptr; + if (dynamic_grid && implicit) { + Jacobian_i = new su2double* [nVar]; + for (auto iVar = 0u; iVar < nVar; iVar++) + Jacobian_i[iVar] = new su2double [nVar] (); + } - /*--- Get coordinates of i & nearest normal and compute distance ---*/ + /*--- Loop over all of the vertices on this boundary marker ---*/ - Coord_i = geometry->nodes->GetCoord(iPoint); - Coord_j = geometry->nodes->GetCoord(Point_Normal); - dist_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) - dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - dist_ij = sqrt(dist_ij); + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - /*--- Store the corrected velocity at the wall which will - be zero (v = 0), unless there is grid motion (v = u_wall)---*/ + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - if (dynamic_grid) { - GridVel = geometry->nodes->GetGridVel(iPoint); - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = GridVel[iDim]; - } - else { - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = 0.0; - } + /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - /*--- Initialize the convective & viscous residuals to zero ---*/ + if (!geometry->nodes->GetDomain(iPoint)) continue; - for (iVar = 0; iVar < nVar; iVar++) { - Res_Conv[iVar] = 0.0; - Res_Visc[iVar] = 0.0; - } + /*--- If it is a customizable patch, retrieve the specified wall heat flux. ---*/ - /*--- Set the residual, truncation error and velocity value on the boundary ---*/ + if (config->GetMarker_All_PyCustom(val_marker)) + Wall_HeatFlux = geometry->GetCustomBoundaryHeatFlux(val_marker, iVertex); - nodes->SetVelocity_Old(iPoint,Vector); + /*--- Compute dual-grid area and boundary normal ---*/ - for (iDim = 0; iDim < nDim; iDim++) - LinSysRes.SetBlock_Zero(iPoint, iDim+1); - nodes->SetVel_ResTruncError_Zero(iPoint); + const auto Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - /*--- Compute the normal gradient in temperature using Twall ---*/ + su2double Area = GeometryToolbox::Norm(nDim, Normal); - dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij; + su2double UnitNormal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + UnitNormal[iDim] = -Normal[iDim]/Area; - /*--- Get transport coefficients ---*/ + /*--- Apply a weak boundary condition for the energy equation. + Compute the residual due to the prescribed heat flux. + The convective part will be zero if the grid is not moving. ---*/ - laminar_viscosity = nodes->GetLaminarViscosity(iPoint); - eddy_viscosity = nodes->GetEddyViscosity(iPoint); - thermal_conductivity = Cp * ( laminar_viscosity/Prandtl_Lam + eddy_viscosity/Prandtl_Turb); + su2double Res_Conv = 0.0; + su2double Res_Visc = Wall_HeatFlux * Area; - // work in progress on real-gases... - //thermal_conductivity = nodes->GetThermalConductivity(iPoint); - //Cp = nodes->GetSpecificHeatCp(iPoint); - //thermal_conductivity += Cp*eddy_viscosity/Prandtl_Turb; + /*--- Impose the value of the velocity as a strong boundary + condition (Dirichlet). Fix the velocity and remove any + contribution to the residual at this node. ---*/ - /*--- Apply a weak boundary condition for the energy equation. - Compute the residual due to the prescribed heat flux. ---*/ + if (dynamic_grid) { + nodes->SetVelocity_Old(iPoint, geometry->nodes->GetGridVel(iPoint)); + } + else { + su2double zero[MAXNDIM] = {0.0}; + nodes->SetVelocity_Old(iPoint, zero); + } - Res_Visc[nDim+1] = thermal_conductivity * dTdn * Area; + for (auto iDim = 0u; iDim < nDim; iDim++) + LinSysRes.SetBlock_Zero(iPoint, iDim+1); + nodes->SetVel_ResTruncError_Zero(iPoint); - /*--- Calculate Jacobian for implicit time stepping ---*/ + /*--- If the wall is moving, there are additional residual contributions + due to pressure (p v_wall.n) and shear stress (tau.v_wall.n). ---*/ + if (dynamic_grid) { if (implicit) { - - for (iVar = 0; iVar < nVar; iVar ++) - for (jVar = 0; jVar < nVar; jVar ++) - Jacobian_i[iVar][jVar] = 0.0; - - /*--- Calculate useful quantities ---*/ - - Density = nodes->GetDensity(iPoint); - Vel2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Vel2 += pow(nodes->GetVelocity(iPoint,iDim),2); - dTdrho = 1.0/Density * ( -Twall + (Gamma-1.0)/Gas_Constant*(Vel2/2.0) ); - - /*--- Enforce the no-slip boundary condition in a strong way ---*/ - - for (iVar = 1; iVar <= nDim; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } - - /*--- Add contributions to the Jacobian from the weak enforcement of the energy equations ---*/ - - Jacobian_i[nDim+1][0] = -thermal_conductivity*theta2/dist_ij * dTdrho * Area; - Jacobian_i[nDim+1][nDim+1] = -thermal_conductivity*theta2/dist_ij * (Gamma-1.0)/(Gas_Constant*Density) * Area; - - /*--- Subtract the block from the Global Jacobian structure ---*/ - - Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); - + for (auto iVar = 0u; iVar < nVar; ++iVar) + Jacobian_i[nDim+1][iVar] = 0.0; } - /*--- If the wall is moving, there are additional residual contributions - due to pressure (p v_wall.n) and shear stress (tau.v_wall.n). ---*/ - - if (dynamic_grid) { - - /*--- Get the grid velocity at the current boundary node ---*/ - - GridVel = geometry->nodes->GetGridVel(iPoint); - ProjGridVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjGridVel += GridVel[iDim]*UnitNormal[iDim]*Area; - - /*--- Retrieve other primitive quantities and viscosities ---*/ - - Density = nodes->GetDensity(iPoint); - Pressure = nodes->GetPressure(iPoint); - laminar_viscosity = nodes->GetLaminarViscosity(iPoint); - eddy_viscosity = nodes->GetEddyViscosity(iPoint); - - total_viscosity = laminar_viscosity + eddy_viscosity; - - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint,iDim+1, jDim); - } - } - - /*--- Divergence of the velocity ---*/ - - div_vel = 0.0; for (iDim = 0 ; iDim < nDim; iDim++) div_vel += Grad_Vel[iDim][iDim]; - - /*--- Compute the viscous stress tensor ---*/ - - for (iDim = 0; iDim < nDim; iDim++) - for (jDim = 0; jDim < nDim; jDim++) { - tau[iDim][jDim] = total_viscosity*( Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim] ) - - TWO3*total_viscosity*div_vel*delta[iDim][jDim]; - } - - /*--- Dot product of the stress tensor with the grid velocity ---*/ - - for (iDim = 0 ; iDim < nDim; iDim++) { - tau_vel[iDim] = 0.0; - for (jDim = 0 ; jDim < nDim; jDim++) - tau_vel[iDim] += tau[iDim][jDim]*GridVel[jDim]; - } - - /*--- Compute the convective and viscous residuals (energy eqn.) ---*/ - - Res_Conv[nDim+1] = Pressure*ProjGridVel; - for (iDim = 0 ; iDim < nDim; iDim++) - Res_Visc[nDim+1] += tau_vel[iDim]*UnitNormal[iDim]*Area; - - /*--- Implicit Jacobian contributions due to moving walls ---*/ + const auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - if (implicit) { - - /*--- Jacobian contribution related to the pressure term ---*/ - - GridVel2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - GridVel2 += GridVel[iDim]*GridVel[iDim]; - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - - Jacobian_i[nDim+1][0] = 0.5*(Gamma-1.0)*GridVel2*ProjGridVel; - for (jDim = 0; jDim < nDim; jDim++) - Jacobian_i[nDim+1][jDim+1] = -(Gamma-1.0)*GridVel[jDim]*ProjGridVel; - Jacobian_i[nDim+1][nDim+1] = (Gamma-1.0)*ProjGridVel; - - /*--- Add the block to the Global Jacobian structure ---*/ - - Jacobian.AddBlock2Diag(iPoint, Jacobian_i); - - /*--- Now the Jacobian contribution related to the shear stress ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - - factor = total_viscosity*Area/(Density*dist_ij); - - if (nDim == 2) { - thetax = theta2 + UnitNormal[0]*UnitNormal[0]/3.0; - thetay = theta2 + UnitNormal[1]*UnitNormal[1]/3.0; - - etaz = UnitNormal[0]*UnitNormal[1]/3.0; - - pix = GridVel[0]*thetax + GridVel[1]*etaz; - piy = GridVel[0]*etaz + GridVel[1]*thetay; + AddDynamicGridResidualContribution(iPoint, Point_Normal, geometry, UnitNormal, + Area, geometry->nodes->GetGridVel(iPoint), + Jacobian_i, Res_Conv, Res_Visc); + } - Jacobian_i[nDim+1][0] -= factor*(-pix*GridVel[0]+piy*GridVel[1]); - Jacobian_i[nDim+1][1] -= factor*pix; - Jacobian_i[nDim+1][2] -= factor*piy; - } - else { - thetax = theta2 + UnitNormal[0]*UnitNormal[0]/3.0; - thetay = theta2 + UnitNormal[1]*UnitNormal[1]/3.0; - thetaz = theta2 + UnitNormal[2]*UnitNormal[2]/3.0; - - etaz = UnitNormal[0]*UnitNormal[1]/3.0; - etax = UnitNormal[1]*UnitNormal[2]/3.0; - etay = UnitNormal[0]*UnitNormal[2]/3.0; - - pix = GridVel[0]*thetax + GridVel[1]*etaz + GridVel[2]*etay; - piy = GridVel[0]*etaz + GridVel[1]*thetay + GridVel[2]*etax; - piz = GridVel[0]*etay + GridVel[1]*etax + GridVel[2]*thetaz; - - Jacobian_i[nDim+1][0] -= factor*(-pix*GridVel[0]+piy*GridVel[1]+piz*GridVel[2]); - Jacobian_i[nDim+1][1] -= factor*pix; - Jacobian_i[nDim+1][2] -= factor*piy; - Jacobian_i[nDim+1][3] -= factor*piz; - } + /*--- Convective and viscous contributions to the residual at the wall ---*/ - /*--- Subtract the block from the Global Jacobian structure ---*/ + LinSysRes(iPoint, nDim+1) += Res_Conv - Res_Visc; - Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); - } + /*--- Enforce the no-slip boundary condition in a strong way by + modifying the velocity-rows of the Jacobian (1 on the diagonal). + And add the contributions to the Jacobian due to energy. ---*/ + if (implicit) { + if (dynamic_grid) { + Jacobian.AddBlock2Diag(iPoint, Jacobian_i); } - /*--- Convective contribution to the residual at the wall ---*/ - - LinSysRes.AddBlock(iPoint, Res_Conv); - - /*--- Viscous contribution to the residual at the wall ---*/ - - LinSysRes.SubtractBlock(iPoint, Res_Visc); - - /*--- Enforce the no-slip boundary condition in a strong way by - modifying the velocity-rows of the Jacobian (1 on the diagonal). ---*/ - - if (implicit) { - for (iVar = 1; iVar <= nDim; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } + for (auto iVar = 1u; iVar <= nDim; iVar++) { + auto total_index = iPoint*nVar+iVar; + Jacobian.DeleteValsRowi(total_index); } - } } -} -void CNSSolver::SetRoe_Dissipation(CGeometry *geometry, CConfig *config){ + if (Jacobian_i) + for (auto iVar = 0u; iVar < nVar; iVar++) + delete [] Jacobian_i[iVar]; + delete [] Jacobian_i; - const unsigned short kind_roe_dissipation = config->GetKind_RoeLowDiss(); +} - SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { +su2double CNSSolver::GetCHTWallTemperature(const CConfig* config, unsigned short val_marker, + unsigned long iVertex, su2double thermal_conductivity, + su2double dist_ij, su2double There, + su2double Temperature_Ref) const { - if (kind_roe_dissipation == FD || kind_roe_dissipation == FD_DUCROS){ + /*--- Compute the normal gradient in temperature using Twall ---*/ - su2double wall_distance = geometry->nodes->GetWall_Distance(iPoint); + const su2double Tconjugate = GetConjugateHeatVariable(val_marker, iVertex, 0) / Temperature_Ref; - nodes->SetRoe_Dissipation_FD(iPoint, wall_distance); + su2double Twall = 0.0; - } else if (kind_roe_dissipation == NTS || kind_roe_dissipation == NTS_DUCROS) { + if ((config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_NEUMANN_HEATFLUX) || + (config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) { - const su2double delta = geometry->nodes->GetMaxLength(iPoint); - assert(delta > 0 && "Delta must be initialized and non-negative"); - nodes->SetRoe_Dissipation_NTS(iPoint, delta, config->GetConst_DES()); - } - } + /*--- Compute wall temperature from both temperatures ---*/ -} + su2double HF_FactorHere = thermal_conductivity*config->GetViscosity_Ref()/dist_ij; + su2double HF_FactorConjugate = GetConjugateHeatVariable(val_marker, iVertex, 2); -void CNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CConfig *config, unsigned short val_marker) { + Twall = (There*HF_FactorHere + Tconjugate*HF_FactorConjugate)/(HF_FactorHere + HF_FactorConjugate); + } + else if ((config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_NEUMANN_HEATFLUX) || + (config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX)) { - unsigned short iVar, jVar, iDim, jDim; - unsigned long iVertex, iPoint, Point_Normal, total_index; + /*--- (Directly) Set wall temperature to conjugate temperature. ---*/ - su2double *Normal, *Coord_i, *Coord_j, Area, dist_ij, theta2; - su2double Twall= 0.0, There, dTdn= 0.0, dTdrho, thermal_conductivity, Tconjugate, HF_FactorHere, HF_FactorConjugate; - su2double thetax, thetay, thetaz, etax, etay, etaz, pix, piy, piz, factor; - su2double ProjGridVel, *GridVel, GridVel2, Pressure = 0.0, Density, Vel2; - su2double total_viscosity, div_vel, tau_vel[3] = {0.0,0.0,0.0}, UnitNormal[3] = {0.0,0.0,0.0}; - su2double laminar_viscosity, eddy_viscosity, Grad_Vel[3][3] = {{1.0, 0.0, 0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}}, - tau[3][3] = {{0.0, 0.0, 0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}}, delta[3][3] = {{1.0, 0.0, 0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}}; + Twall = Tconjugate; + } + else { + SU2_MPI::Error("Unknown CHT coupling method.", CURRENT_FUNCTION); + } - su2double Prandtl_Lam = config->GetPrandtl_Lam(); - su2double Prandtl_Turb = config->GetPrandtl_Turb(); - su2double Gas_Constant = config->GetGas_ConstantND(); - su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; + return Twall; +} - su2double Temperature_Ref = config->GetTemperature_Ref(); +void CNSSolver::BC_Isothermal_Wall_Generic(CGeometry *geometry, CSolver **solver_container, + CNumerics *conv_numerics, CNumerics *visc_numerics, + CConfig *config, unsigned short val_marker, bool cht_mode) { - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const su2double Temperature_Ref = config->GetTemperature_Ref(); + const su2double Prandtl_Lam = config->GetPrandtl_Lam(); + const su2double Prandtl_Turb = config->GetPrandtl_Turb(); + const su2double Gas_Constant = config->GetGas_ConstantND(); + const su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; - /*--- Identify the boundary ---*/ + /*--- Identify the boundary and retrieve the specified wall temperature from + the config (for non-CHT problems) as well as the wall function treatment. ---*/ - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); + su2double Twall = 0.0; + if (!cht_mode) { + Twall = config->GetIsothermal_Temperature(Marker_Tag) / Temperature_Ref; + } -// /*--- Retrieve the specified wall function treatment.---*/ -// // Wall_Function = config->GetWallFunction_Treatment(Marker_Tag); // if (Wall_Function != NO_WALL_FUNCTION) { -// SU2_MPI::Error("Wall function treament not implemented yet", CURRENT_FUNCTION); +// SU2_MPI::Error("Wall function treament not implemented yet", CURRENT_FUNCTION); // } - /*--- Loop over boundary points ---*/ - - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Compute dual-grid area and boundary normal ---*/ - - Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - - Area = 0.0; for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; Area = sqrt (Area); - - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = -Normal[iDim]/Area; - - /*--- Calculate useful quantities ---*/ - - theta2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - theta2 += UnitNormal[iDim]*UnitNormal[iDim]; - - /*--- Compute closest normal neighbor ---*/ - - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - - /*--- Get coordinates of i & nearest normal and compute distance ---*/ - - Coord_i = geometry->nodes->GetCoord(iPoint); - Coord_j = geometry->nodes->GetCoord(Point_Normal); - dist_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) - dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - dist_ij = sqrt(dist_ij); - - /*--- Store the corrected velocity at the wall which will - be zero (v = 0), unless there is grid motion (v = u_wall)---*/ - - if (dynamic_grid) { - GridVel = geometry->nodes->GetGridVel(iPoint); - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = GridVel[iDim]; - } - else { - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = 0.0; - } - - /*--- Initialize the convective & viscous residuals to zero ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Res_Conv[iVar] = 0.0; - Res_Visc[iVar] = 0.0; - } - - /*--- Set the residual, truncation error and velocity value on the boundary ---*/ - - nodes->SetVelocity_Old(iPoint,Vector); - - for (iDim = 0; iDim < nDim; iDim++) - LinSysRes.SetBlock_Zero(iPoint, iDim+1); - nodes->SetVel_ResTruncError_Zero(iPoint); - - /*--- Get transport coefficients ---*/ - - laminar_viscosity = nodes->GetLaminarViscosity(iPoint); - eddy_viscosity = nodes->GetEddyViscosity(iPoint); - thermal_conductivity = Cp * ( laminar_viscosity/Prandtl_Lam + eddy_viscosity/Prandtl_Turb); - - // work in progress on real-gases... - //thermal_conductivity = nodes->GetThermalConductivity(iPoint); - //Cp = nodes->GetSpecificHeatCp(iPoint); - //thermal_conductivity += Cp*eddy_viscosity/Prandtl_Turb; - - /*--- Compute the normal gradient in temperature using Twall ---*/ - - There = nodes->GetTemperature(Point_Normal); - Tconjugate = GetConjugateHeatVariable(val_marker, iVertex, 0)/Temperature_Ref; - - if ((config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_NEUMANN_HEATFLUX) || - (config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) { - - /*--- Compute wall temperature from both temperatures ---*/ - - HF_FactorHere = thermal_conductivity*config->GetViscosity_Ref()/dist_ij; - HF_FactorConjugate = GetConjugateHeatVariable(val_marker, iVertex, 2); - - Twall = (There*HF_FactorHere + Tconjugate*HF_FactorConjugate)/(HF_FactorHere + HF_FactorConjugate); - dTdn = -(There - Twall)/dist_ij; - } - else if ((config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_NEUMANN_HEATFLUX) || - (config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX)) { - - /*--- (Directly) Set wall temperature to conjugate temperature. ---*/ - - Twall = Tconjugate; - dTdn = -(There - Twall)/dist_ij; - } - else { - Twall = dTdn = 0.0; - SU2_MPI::Error("Unknown CHT coupling method.", CURRENT_FUNCTION); - } - - /*--- Apply a weak boundary condition for the energy equation. - Compute the residual due to the prescribed heat flux. ---*/ - - Res_Visc[nDim+1] = thermal_conductivity * dTdn * Area; - - /*--- Calculate Jacobian for implicit time stepping ---*/ + su2double **Jacobian_i = nullptr; + if (implicit) { + Jacobian_i = new su2double* [nVar]; + for (auto iVar = 0u; iVar < nVar; iVar++) + Jacobian_i[iVar] = new su2double [nVar] (); + } - if (implicit) { + /*--- Loop over boundary points ---*/ - for (iVar = 0; iVar < nVar; iVar ++) - for (jVar = 0; jVar < nVar; jVar ++) - Jacobian_i[iVar][jVar] = 0.0; + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - /*--- Calculate useful quantities ---*/ + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - Density = nodes->GetDensity(iPoint); - Vel2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Vel2 += pow(nodes->GetVelocity(iPoint,iDim),2); - dTdrho = 1.0/Density * ( -Twall + (Gamma-1.0)/Gas_Constant*(Vel2/2.0) ); + if (!geometry->nodes->GetDomain(iPoint)) continue; - /*--- Enforce the no-slip boundary condition in a strong way ---*/ + /*--- Compute dual-grid area and boundary normal ---*/ - for (iVar = 1; iVar <= nDim; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } + const auto Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - /*--- Add contributions to the Jacobian from the weak enforcement of the energy equations ---*/ + su2double Area = GeometryToolbox::Norm(nDim, Normal); - Jacobian_i[nDim+1][0] = -thermal_conductivity*theta2/dist_ij * dTdrho * Area; - Jacobian_i[nDim+1][nDim+1] = -thermal_conductivity*theta2/dist_ij * (Gamma-1.0)/(Gas_Constant*Density) * Area; + su2double UnitNormal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + UnitNormal[iDim] = -Normal[iDim]/Area; - /*--- Subtract the block from the Global Jacobian structure ---*/ + /*--- Compute closest normal neighbor ---*/ - Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); + const auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - } + /*--- Get coordinates of i & nearest normal and compute distance ---*/ - /*--- If the wall is moving, there are additional residual contributions - due to pressure (p v_wall.n) and shear stress (tau.v_wall.n). ---*/ + const auto Coord_i = geometry->nodes->GetCoord(iPoint); + const auto Coord_j = geometry->nodes->GetCoord(Point_Normal); - if (dynamic_grid) { + su2double dist_ij = GeometryToolbox::Distance(nDim, Coord_i, Coord_j); - /*--- Get the grid velocity at the current boundary node ---*/ + /*--- Store the corrected velocity at the wall which will + be zero (v = 0), unless there is grid motion (v = u_wall)---*/ - GridVel = geometry->nodes->GetGridVel(iPoint); - ProjGridVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjGridVel += GridVel[iDim]*UnitNormal[iDim]*Area; - - /*--- Retrieve other primitive quantities and viscosities ---*/ - - Density = nodes->GetDensity(iPoint); - Pressure = nodes->GetPressure(iPoint); - laminar_viscosity = nodes->GetLaminarViscosity(iPoint); - eddy_viscosity = nodes->GetEddyViscosity(iPoint); - - total_viscosity = laminar_viscosity + eddy_viscosity; - - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - Grad_Vel[iDim][jDim] = nodes->GetGradient_Primitive(iPoint,iDim+1, jDim); - } - } + if (dynamic_grid) { + nodes->SetVelocity_Old(iPoint, geometry->nodes->GetGridVel(iPoint)); + } + else { + su2double zero[MAXNDIM] = {0.0}; + nodes->SetVelocity_Old(iPoint, zero); + } - /*--- Divergence of the velocity ---*/ + for (auto iDim = 0u; iDim < nDim; iDim++) + LinSysRes.SetBlock_Zero(iPoint, iDim+1); + nodes->SetVel_ResTruncError_Zero(iPoint); - div_vel = 0.0; for (iDim = 0 ; iDim < nDim; iDim++) div_vel += Grad_Vel[iDim][iDim]; + /*--- Get transport coefficients ---*/ - /*--- Compute the viscous stress tensor ---*/ + su2double laminar_viscosity = nodes->GetLaminarViscosity(iPoint); + su2double eddy_viscosity = nodes->GetEddyViscosity(iPoint); + su2double thermal_conductivity = Cp * (laminar_viscosity/Prandtl_Lam + eddy_viscosity/Prandtl_Turb); - for (iDim = 0; iDim < nDim; iDim++) - for (jDim = 0; jDim < nDim; jDim++) { - tau[iDim][jDim] = total_viscosity*( Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim] ) - - TWO3*total_viscosity*div_vel*delta[iDim][jDim]; - } + // work in progress on real-gases... + //thermal_conductivity = nodes->GetThermalConductivity(iPoint); + //Cp = nodes->GetSpecificHeatCp(iPoint); + //thermal_conductivity += Cp*eddy_viscosity/Prandtl_Turb; - /*--- Dot product of the stress tensor with the grid velocity ---*/ + /*--- If it is a customizable or CHT patch, retrieve the specified wall temperature. ---*/ - for (iDim = 0 ; iDim < nDim; iDim++) { - tau_vel[iDim] = 0.0; - for (jDim = 0 ; jDim < nDim; jDim++) - tau_vel[iDim] += tau[iDim][jDim]*GridVel[jDim]; - } + const su2double There = nodes->GetTemperature(Point_Normal); - /*--- Compute the convective and viscous residuals (energy eqn.) ---*/ + if (cht_mode) { + Twall = GetCHTWallTemperature(config, val_marker, iVertex, dist_ij, + thermal_conductivity, There, Temperature_Ref); + } + else if (config->GetMarker_All_PyCustom(val_marker)) { + Twall = geometry->GetCustomBoundaryTemperature(val_marker, iVertex); + } - Res_Conv[nDim+1] = Pressure*ProjGridVel; - for (iDim = 0 ; iDim < nDim; iDim++) - Res_Visc[nDim+1] += tau_vel[iDim]*UnitNormal[iDim]*Area; + /*--- Compute the normal gradient in temperature using Twall ---*/ - /*--- Implicit Jacobian contributions due to moving walls ---*/ + su2double dTdn = -(There - Twall)/dist_ij; - if (implicit) { + /*--- Apply a weak boundary condition for the energy equation. + Compute the residual due to the prescribed heat flux. ---*/ - /*--- Jacobian contribution related to the pressure term ---*/ + su2double Res_Conv = 0.0; + su2double Res_Visc = thermal_conductivity * dTdn * Area; - GridVel2 = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - GridVel2 += GridVel[iDim]*GridVel[iDim]; - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; + /*--- Calculate Jacobian for implicit time stepping ---*/ - Jacobian_i[nDim+1][0] = 0.5*(Gamma-1.0)*GridVel2*ProjGridVel; - for (jDim = 0; jDim < nDim; jDim++) - Jacobian_i[nDim+1][jDim+1] = -(Gamma-1.0)*GridVel[jDim]*ProjGridVel; - Jacobian_i[nDim+1][nDim+1] = (Gamma-1.0)*ProjGridVel; + if (implicit) { - /*--- Add the block to the Global Jacobian structure ---*/ + /*--- Add contributions to the Jacobian from the weak enforcement of the energy equations. ---*/ - Jacobian.AddBlock2Diag(iPoint, Jacobian_i); + su2double Density = nodes->GetDensity(iPoint); + su2double Vel2 = GeometryToolbox::SquaredNorm(nDim, &nodes->GetPrimitive(iPoint)[1]); + su2double dTdrho = 1.0/Density * ( -Twall + (Gamma-1.0)/Gas_Constant*(Vel2/2.0) ); - /*--- Now the Jacobian contribution related to the shear stress ---*/ + Jacobian_i[nDim+1][0] = thermal_conductivity/dist_ij * dTdrho * Area; - for (iVar = 0; iVar < nVar; iVar++) - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; + for (auto jDim = 0u; jDim < nDim; jDim++) + Jacobian_i[nDim+1][jDim+1] = 0.0; - factor = total_viscosity*Area/(Density*dist_ij); + Jacobian_i[nDim+1][nDim+1] = thermal_conductivity/dist_ij * (Gamma-1.0)/(Gas_Constant*Density) * Area; + } - if (nDim == 2) { - thetax = theta2 + UnitNormal[0]*UnitNormal[0]/3.0; - thetay = theta2 + UnitNormal[1]*UnitNormal[1]/3.0; + /*--- If the wall is moving, there are additional residual contributions + due to pressure (p v_wall.n) and shear stress (tau.v_wall.n). ---*/ - etaz = UnitNormal[0]*UnitNormal[1]/3.0; + if (dynamic_grid) { + AddDynamicGridResidualContribution(iPoint, Point_Normal, geometry, UnitNormal, + Area, geometry->nodes->GetGridVel(iPoint), + Jacobian_i, Res_Conv, Res_Visc); + } - pix = GridVel[0]*thetax + GridVel[1]*etaz; - piy = GridVel[0]*etaz + GridVel[1]*thetay; + /*--- Convective and viscous contributions to the residual at the wall ---*/ - Jacobian_i[nDim+1][0] -= factor*(-pix*GridVel[0]+piy*GridVel[1]); - Jacobian_i[nDim+1][1] -= factor*pix; - Jacobian_i[nDim+1][2] -= factor*piy; - } - else { - thetax = theta2 + UnitNormal[0]*UnitNormal[0]/3.0; - thetay = theta2 + UnitNormal[1]*UnitNormal[1]/3.0; - thetaz = theta2 + UnitNormal[2]*UnitNormal[2]/3.0; - - etaz = UnitNormal[0]*UnitNormal[1]/3.0; - etax = UnitNormal[1]*UnitNormal[2]/3.0; - etay = UnitNormal[0]*UnitNormal[2]/3.0; - - pix = GridVel[0]*thetax + GridVel[1]*etaz + GridVel[2]*etay; - piy = GridVel[0]*etaz + GridVel[1]*thetay + GridVel[2]*etax; - piz = GridVel[0]*etay + GridVel[1]*etax + GridVel[2]*thetaz; - - Jacobian_i[nDim+1][0] -= factor*(-pix*GridVel[0]+piy*GridVel[1]+piz*GridVel[2]); - Jacobian_i[nDim+1][1] -= factor*pix; - Jacobian_i[nDim+1][2] -= factor*piy; - Jacobian_i[nDim+1][3] -= factor*piz; - } + LinSysRes(iPoint, nDim+1) += Res_Conv - Res_Visc; - /*--- Subtract the block from the Global Jacobian structure ---*/ + /*--- Enforce the no-slip boundary condition in a strong way by + modifying the velocity-rows of the Jacobian (1 on the diagonal). + And add the contributions to the Jacobian due to energy. ---*/ - Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); - } + if (implicit) { + Jacobian.AddBlock2Diag(iPoint, Jacobian_i); + for (auto iVar = 1u; iVar <= nDim; iVar++) { + auto total_index = iPoint*nVar+iVar; + Jacobian.DeleteValsRowi(total_index); } + } + } - /*--- Convective contribution to the residual at the wall ---*/ + if (Jacobian_i) + for (auto iVar = 0u; iVar < nVar; iVar++) + delete [] Jacobian_i[iVar]; + delete [] Jacobian_i; - LinSysRes.AddBlock(iPoint, Res_Conv); +} - /*--- Viscous contribution to the residual at the wall ---*/ +void CNSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - LinSysRes.SubtractBlock(iPoint, Res_Visc); + BC_Isothermal_Wall_Generic(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); +} - /*--- Enforce the no-slip boundary condition in a strong way by - modifying the velocity-rows of the Jacobian (1 on the diagonal). ---*/ +void CNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CConfig *config, unsigned short val_marker) { - if (implicit) { - for (iVar = 1; iVar <= nDim; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } - } - } - } + BC_Isothermal_Wall_Generic(geometry, solver_container, conv_numerics, nullptr, config, val_marker, true); } void CNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container, CConfig *config) { - unsigned short iDim, jDim, iMarker; - unsigned long iVertex, iPoint, Point_Normal, counter; - - su2double Area, div_vel, UnitNormal[3], *Normal; - su2double **grad_primvar, tau[3][3]; + const su2double Gas_Constant = config->GetGas_ConstantND(); + const su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; - su2double Vel[3] = {0.0, 0.0, 0.0}, VelNormal, VelTang[3], VelTangMod, VelInfMod, WallDist[3], WallDistMod; - su2double T_Normal, P_Normal; - su2double Density_Wall, T_Wall, P_Wall, Lam_Visc_Wall, Tau_Wall = 0.0, Tau_Wall_Old = 0.0; - su2double *Coord, *Coord_Normal; - su2double diff, Delta; - su2double U_Tau, U_Plus, Gam, Beta, Phi, Q, Y_Plus_White, Y_Plus; - su2double TauElem[3], TauNormal, TauTangent[3], WallShearStress; - su2double Gas_Constant = config->GetGas_ConstantND(); - su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; - - unsigned short max_iter = 10; - su2double tol = 1e-6; - - /*--- Get the freestream velocity magnitude for non-dim. purposes ---*/ - - su2double *VelInf = config->GetVelocity_FreeStreamND(); - VelInfMod = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - VelInfMod += VelInf[iDim]; - VelInfMod = sqrt(VelInfMod); + constexpr unsigned short max_iter = 10; + const su2double tol = 1e-6; + const su2double relax = 0.25; /*--- Compute the recovery factor ---*/ // Double-check: laminar or turbulent Pr for this? - su2double Recovery = pow(config->GetPrandtl_Lam(), (1.0/3.0)); + const su2double Recovery = pow(config->GetPrandtl_Lam(), (1.0/3.0)); /*--- Typical constants from boundary layer theory ---*/ - su2double kappa = 0.4; - su2double B = 5.5; + const su2double kappa = 0.4; + const su2double B = 5.5; - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) == HEAT_FLUX) || - (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) ) { + if (!config->GetViscous_Wall(iMarker)) continue; - /*--- Identify the boundary by string name ---*/ + /*--- Identify the boundary by string name ---*/ - string Marker_Tag = config->GetMarker_All_TagBound(iMarker); + const auto Marker_Tag = config->GetMarker_All_TagBound(iMarker); - /*--- Get the specified wall heat flux from config ---*/ + /*--- Get the specified wall heat flux from config ---*/ - // Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag); + // Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag); - /*--- Loop over all of the vertices on this boundary marker ---*/ + /*--- Loop over all of the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + SU2_OMP_FOR_DYN(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - Point_Normal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor(); - - /*--- Check if the node belongs to the domain (i.e, not a halo node) - and the neighbor is not part of the physical boundary ---*/ - - if (geometry->nodes->GetDomain(iPoint)) { + const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + const auto Point_Normal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor(); - /*--- Get coordinates of the current vertex and nearest normal point ---*/ + /*--- Check if the node belongs to the domain (i.e, not a halo node) + and the neighbor is not part of the physical boundary ---*/ - Coord = geometry->nodes->GetCoord(iPoint); - Coord_Normal = geometry->nodes->GetCoord(Point_Normal); + if (!geometry->nodes->GetDomain(iPoint)) continue; - /*--- Compute dual-grid area and boundary normal ---*/ - - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + /*--- Get coordinates of the current vertex and nearest normal point ---*/ - Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Area += Normal[iDim]*Normal[iDim]; - Area = sqrt (Area); + const auto Coord = geometry->nodes->GetCoord(iPoint); + const auto Coord_Normal = geometry->nodes->GetCoord(Point_Normal); - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = -Normal[iDim]/Area; + /*--- Compute dual-grid area and boundary normal ---*/ - /*--- Get the velocity, pressure, and temperature at the nearest - (normal) interior point. ---*/ + const auto Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - for (iDim = 0; iDim < nDim; iDim++) - Vel[iDim] = nodes->GetVelocity(Point_Normal,iDim); - P_Normal = nodes->GetPressure(Point_Normal); - T_Normal = nodes->GetTemperature(Point_Normal); + su2double Area = GeometryToolbox::Norm(nDim, Normal); - /*--- Compute the wall-parallel velocity at first point off the wall ---*/ + su2double UnitNormal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + UnitNormal[iDim] = -Normal[iDim]/Area; - VelNormal = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - VelNormal += Vel[iDim] * UnitNormal[iDim]; - for (iDim = 0; iDim < nDim; iDim++) - VelTang[iDim] = Vel[iDim] - VelNormal*UnitNormal[iDim]; + /*--- Get the velocity, pressure, and temperature at the nearest + (normal) interior point. ---*/ - VelTangMod = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - VelTangMod += VelTang[iDim]*VelTang[iDim]; - VelTangMod = sqrt(VelTangMod); + su2double Vel[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Vel[iDim] = nodes->GetVelocity(Point_Normal,iDim); + su2double P_Normal = nodes->GetPressure(Point_Normal); + su2double T_Normal = nodes->GetTemperature(Point_Normal); - /*--- Compute normal distance of the interior point from the wall ---*/ + /*--- Compute the wall-parallel velocity at first point off the wall ---*/ - for (iDim = 0; iDim < nDim; iDim++) - WallDist[iDim] = (Coord[iDim] - Coord_Normal[iDim]); + su2double VelNormal = GeometryToolbox::DotProduct(nDim, Vel, UnitNormal); - WallDistMod = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - WallDistMod += WallDist[iDim]*WallDist[iDim]; - WallDistMod = sqrt(WallDistMod); + su2double VelTang[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + VelTang[iDim] = Vel[iDim] - VelNormal*UnitNormal[iDim]; - /*--- Compute mach number ---*/ + su2double VelTangMod = GeometryToolbox::Norm(int(MAXNDIM), VelTang); - // M_Normal = VelTangMod / sqrt(Gamma * Gas_Constant * T_Normal); + /*--- Compute normal distance of the interior point from the wall ---*/ - /*--- Compute the wall temperature using the Crocco-Buseman equation ---*/ + su2double WallDist[MAXNDIM] = {0.0}; + GeometryToolbox::Distance(nDim, Coord, Coord_Normal, WallDist); - //T_Wall = T_Normal * (1.0 + 0.5*Gamma_Minus_One*Recovery*M_Normal*M_Normal); - T_Wall = T_Normal + Recovery*pow(VelTangMod,2.0)/(2.0*Cp); + su2double WallDistMod = GeometryToolbox::Norm(int(MAXNDIM), WallDist); - /*--- Extrapolate the pressure from the interior & compute the - wall density using the equation of state ---*/ + /*--- Compute mach number ---*/ - P_Wall = P_Normal; - Density_Wall = P_Wall/(Gas_Constant*T_Wall); + // M_Normal = VelTangMod / sqrt(Gamma * Gas_Constant * T_Normal); - /*--- Compute the shear stress at the wall in the regular fashion - by using the stress tensor on the surface ---*/ + /*--- Compute the wall temperature using the Crocco-Buseman equation ---*/ - Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint); - grad_primvar = nodes->GetGradient_Primitive(iPoint); + //T_Wall = T_Normal * (1.0 + 0.5*Gamma_Minus_One*Recovery*M_Normal*M_Normal); + su2double T_Wall = T_Normal + Recovery*pow(VelTangMod,2.0)/(2.0*Cp); - div_vel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - div_vel += grad_primvar[iDim+1][iDim]; + /*--- Extrapolate the pressure from the interior & compute the + wall density using the equation of state ---*/ - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - Delta = 0.0; if (iDim == jDim) Delta = 1.0; - tau[iDim][jDim] = Lam_Visc_Wall*( grad_primvar[jDim+1][iDim] - + grad_primvar[iDim+1][jDim]) - - TWO3*Lam_Visc_Wall*div_vel*Delta; - } - TauElem[iDim] = 0.0; - for (jDim = 0; jDim < nDim; jDim++) - TauElem[iDim] += tau[iDim][jDim]*UnitNormal[jDim]; - } + su2double P_Wall = P_Normal; + su2double Density_Wall = P_Wall/(Gas_Constant*T_Wall); - /*--- Compute wall shear stress as the magnitude of the wall-tangential - component of the shear stress tensor---*/ + /*--- Compute the shear stress at the wall in the regular fashion + by using the stress tensor on the surface ---*/ - TauNormal = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - TauNormal += TauElem[iDim] * UnitNormal[iDim]; + su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint); - for (iDim = 0; iDim < nDim; iDim++) - TauTangent[iDim] = TauElem[iDim] - TauNormal * UnitNormal[iDim]; + const auto GradVel = &nodes->GetGradient_Primitive(iPoint)[1]; - WallShearStress = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - WallShearStress += TauTangent[iDim]*TauTangent[iDim]; - WallShearStress = sqrt(WallShearStress); + su2double div_vel = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) + div_vel += GradVel[iDim][iDim]; - /*--- Calculate the quantities from boundary layer theory and - iteratively solve for a new wall shear stress. Use the current wall - shear stress as a starting guess for the wall function. ---*/ + su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}, TauElem[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) { + for (auto jDim = 0u; jDim < nDim; jDim++) { + tau[iDim][jDim] = Lam_Visc_Wall * (GradVel[jDim][iDim] + GradVel[iDim][jDim]); + } + tau[iDim][iDim] -= TWO3*Lam_Visc_Wall*div_vel; - Tau_Wall_Old = WallShearStress; - counter = 0; diff = 1.0; + TauElem[iDim] = GeometryToolbox::DotProduct(nDim, tau[iDim], UnitNormal); + } - while (diff > tol) { + /*--- Compute wall shear stress as the magnitude of the wall-tangential + component of the shear stress tensor---*/ - /*--- Friction velocity and u+ ---*/ + su2double TauNormal = GeometryToolbox::DotProduct(nDim, TauElem, UnitNormal); - U_Tau = sqrt(Tau_Wall_Old/Density_Wall); - U_Plus = VelTangMod/U_Tau; + su2double TauTangent[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + TauTangent[iDim] = TauElem[iDim] - TauNormal * UnitNormal[iDim]; - /*--- Gamma, Beta, Q, and Phi, defined by Nichols & Nelson (2004) ---*/ + su2double Tau_Wall = GeometryToolbox::Norm(int(MAXNDIM), TauTangent); - Gam = Recovery*U_Tau*U_Tau/(2.0*Cp*T_Wall); - Beta = 0.0; // For adiabatic flows only - Q = sqrt(Beta*Beta + 4.0*Gam); - Phi = asin(-1.0*Beta/Q); + /*--- Calculate the quantities from boundary layer theory and + iteratively solve for a new wall shear stress. Use the current wall + shear stress as a starting guess for the wall function. ---*/ - /*--- Y+ defined by White & Christoph (compressibility and heat transfer) negative value for (2.0*Gam*U_Plus - Beta)/Q ---*/ + su2double Tau_Wall_Old = Tau_Wall; + unsigned short counter = 0; + su2double diff = 1.0; - Y_Plus_White = exp((kappa/sqrt(Gam))*(asin((2.0*Gam*U_Plus - Beta)/Q) - Phi))*exp(-1.0*kappa*B); + while (diff > tol) { - /*--- Spalding's universal form for the BL velocity with the - outer velocity form of White & Christoph above. ---*/ + /*--- Friction velocity and u+ ---*/ - Y_Plus = U_Plus + Y_Plus_White - (exp(-1.0*kappa*B)* - (1.0 + kappa*U_Plus + kappa*kappa*U_Plus*U_Plus/2.0 + - kappa*kappa*kappa*U_Plus*U_Plus*U_Plus/6.0)); + su2double U_Tau = sqrt(Tau_Wall_Old/Density_Wall); + su2double U_Plus = VelTangMod/U_Tau; - /*--- Calculate an updated value for the wall shear stress - using the y+ value, the definition of y+, and the definition of - the friction velocity. ---*/ + /*--- Gamma, Beta, Q, and Phi, defined by Nichols & Nelson (2004) ---*/ - Tau_Wall = (1.0/Density_Wall)*pow(Y_Plus*Lam_Visc_Wall/WallDistMod,2.0); + su2double Gam = Recovery*pow(U_Tau,2)/(2.0*Cp*T_Wall); + su2double Beta = 0.0; // For adiabatic flows only + su2double Q = sqrt(Beta*Beta + 4.0*Gam); + su2double Phi = asin(-1.0*Beta/Q); - /*--- Difference between the old and new Tau. Update old value. ---*/ + /*--- Y+ defined by White & Christoph (compressibility and heat transfer) + negative value for (2.0*Gam*U_Plus - Beta)/Q ---*/ - diff = fabs(Tau_Wall-Tau_Wall_Old); - Tau_Wall_Old += 0.25*(Tau_Wall-Tau_Wall_Old); + su2double Y_Plus_White = exp((kappa/sqrt(Gam))*(asin((2.0*Gam*U_Plus - Beta)/Q) - Phi))*exp(-1.0*kappa*B); - counter++; - if (counter > max_iter) { - cout << "WARNING: Tau_Wall evaluation has not converged in solver_direct_mean.cpp" << endl; - cout << Tau_Wall_Old << " " << Tau_Wall << " " << diff << endl; - break; - } + /*--- Spalding's universal form for the BL velocity with the + outer velocity form of White & Christoph above. ---*/ - } + su2double kUp = kappa*U_Plus; + su2double Y_Plus = U_Plus + Y_Plus_White - exp(-1.0*kappa*B) * (1.0 + kUp*(1.0 + 0.5*kUp + pow(kUp,2)/6.0)); + /*--- Calculate an updated value for the wall shear stress using the y+ value, + the definition of y+, and the definition of the friction velocity. ---*/ - /*--- Store this value for the wall shear stress at the node. ---*/ + Tau_Wall = (1.0/Density_Wall)*pow(Y_Plus*Lam_Visc_Wall/WallDistMod,2.0); - nodes->SetTauWall(iPoint,Tau_Wall); + /*--- Difference between the old and new Tau. Update old value. ---*/ + diff = fabs(Tau_Wall-Tau_Wall_Old); + Tau_Wall_Old += relax * (Tau_Wall-Tau_Wall_Old); + counter++; + if (counter > max_iter) { + cout << "WARNING: Tau_Wall evaluation has not converged in CNSSolver.cpp" << endl; + cout << Tau_Wall_Old << " " << Tau_Wall << " " << diff << endl; + break; } - } + /*--- Store this value for the wall shear stress at the node. ---*/ + + nodes->SetTauWall(iPoint, Tau_Wall); + } + } } diff --git a/SU2_CFD/src/solvers/CRadP1Solver.cpp b/SU2_CFD/src/solvers/CRadP1Solver.cpp index 33bf8b0bd451..d6c132f6692b 100644 --- a/SU2_CFD/src/solvers/CRadP1Solver.cpp +++ b/SU2_CFD/src/solvers/CRadP1Solver.cpp @@ -120,21 +120,6 @@ CRadP1Solver::CRadP1Solver(CGeometry* geometry, CConfig *config) : CRadSolver(ge } } - /*--- Define some auxiliary vectors for computing flow variable - gradients by least squares, S matrix := inv(R)*traspose(inv(R)), - c vector := transpose(WA)*(Wb) ---*/ - - if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { - - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - /*--- c vector := transpose(WA)*(Wb) ---*/ - Cvector = new su2double* [nVar+1]; - for (iVar = 0; iVar < nVar+1; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - /*--- Always instantiate and initialize the variable to a zero value. ---*/ su2double init_val; switch(config->GetKind_P1_Init()){ diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 3ccd3f4d37f7..079ac59aaa50 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -44,6 +44,7 @@ #include "../../../Common/include/toolboxes/MMS/CUserDefinedSolution.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" #include "../../../Common/include/toolboxes/C1DInterpolation.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../include/CMarkerProfileReaderFVM.hpp" @@ -93,8 +94,6 @@ CSolver::CSolver(bool mesh_deform_mode) : System(mesh_deform_mode) { Jacobian_jj = nullptr; iPoint_UndLapl = nullptr; jPoint_UndLapl = nullptr; - Smatrix = nullptr; - Cvector = nullptr; Restart_Vars = nullptr; Restart_Data = nullptr; base_nodes = nullptr; @@ -139,14 +138,12 @@ CSolver::CSolver(bool mesh_deform_mode) : System(mesh_deform_mode) { CSolver::~CSolver(void) { - unsigned short iVar, iDim; + unsigned short iVar; unsigned long iMarker, iVertex; /*--- Public variables, may be accessible outside ---*/ - - delete [] OutputHeadingNames; - + delete [] OutputHeadingNames; /*--- Private ---*/ @@ -227,18 +224,6 @@ CSolver::~CSolver(void) { delete [] Jacobian_jj; } - if (Smatrix != nullptr) { - for (iDim = 0; iDim < nDim; iDim++) - delete [] Smatrix[iDim]; - delete [] Smatrix; - } - - if (Cvector != nullptr) { - for (iVar = 0; iVar < nVarGrad; iVar++) - delete [] Cvector[iVar]; - delete [] Cvector; - } - if (VertexTraction != nullptr) { for (iMarker = 0; iMarker < nMarker; iMarker++) { for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) @@ -259,10 +244,10 @@ CSolver::~CSolver(void) { delete [] nVertex; - if (Restart_Vars != nullptr) {delete [] Restart_Vars; Restart_Vars = nullptr;} - if (Restart_Data != nullptr) {delete [] Restart_Data; Restart_Data = nullptr;} + delete [] Restart_Vars; + delete [] Restart_Data; - if (VerificationSolution != nullptr) {delete VerificationSolution; VerificationSolution = nullptr;} + delete VerificationSolution; } @@ -304,6 +289,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, su2double rotMatrix[3][3] = {{1.0,0.0,0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}}; su2double Theta, Phi, Psi, cosTheta, sinTheta, cosPhi, sinPhi, cosPsi, sinPsi; su2double rotCoord_i[3] = {0.0, 0.0, 0.0}, rotCoord_j[3] = {0.0, 0.0, 0.0}; + su2activematrix Cvector; string Marker_Tag; @@ -934,9 +920,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Inizialization of variables ---*/ - for (iVar = 0; iVar < nVar; iVar++) - for (iDim = 0; iDim < nDim; iDim++) - Cvector[iVar][iDim] = 0.0; + Cvector.resize(nVar,nDim) = su2double(0.0); r11 = 0.0; r12 = 0.0; r22 = 0.0; r13 = 0.0; r23_a = 0.0; r23_b = 0.0; r33 = 0.0; @@ -1035,7 +1019,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) for (iDim = 0; iDim < nDim; iDim++) - Cvector[iVar][iDim] += ((rotCoord_j[iDim]-rotCoord_i[iDim])* + Cvector(iVar,iDim) += ((rotCoord_j[iDim]-rotCoord_i[iDim])* (rotPrim_j[iVar]-rotPrim_i[iVar])/weight); } @@ -1069,7 +1053,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, for (iVar = 0; iVar < nVar; iVar++) { for (iDim = 0; iDim < nDim; iDim++) { - bufDSend[buf_offset] = Cvector[iVar][iDim]; + bufDSend[buf_offset] = Cvector(iVar,iDim); buf_offset++; } } @@ -1143,9 +1127,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Inizialization of variables ---*/ - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - for (iDim = 0; iDim < nDim; iDim++) - Cvector[iVar][iDim] = 0.0; + Cvector.resize(nPrimVarGrad,nDim) = su2double(0.0); r11 = 0.0; r12 = 0.0; r22 = 0.0; r13 = 0.0; r23_a = 0.0; r23_b = 0.0; r33 = 0.0; @@ -1244,7 +1226,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, for (iVar = 0; iVar < nPrimVarGrad; iVar++) for (iDim = 0; iDim < nDim; iDim++) - Cvector[iVar][iDim] += ((rotCoord_j[iDim]-rotCoord_i[iDim])* + Cvector(iVar,iDim) += ((rotCoord_j[iDim]-rotCoord_i[iDim])* (rotPrim_j[iVar]-rotPrim_i[iVar])/weight); } @@ -1278,7 +1260,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, for (iVar = 0; iVar < nPrimVarGrad; iVar++) { for (iDim = 0; iDim < nDim; iDim++) { - bufDSend[buf_offset] = Cvector[iVar][iDim]; + bufDSend[buf_offset] = Cvector(iVar,iDim); buf_offset++; } } @@ -1538,6 +1520,13 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, su2double Time_Step, Volume, Solution_Min, Solution_Max, Limiter_Min; + su2double **Jacobian_i = nullptr; + if ((commType == PERIODIC_RESIDUAL) && implicit_periodic) { + Jacobian_i = new su2double* [nVar]; + for (iVar = 0; iVar < nVar; iVar++) + Jacobian_i[iVar] = new su2double [nVar]; + } + /*--- Set some local pointers to make access simpler. ---*/ su2double *bufDRecv = geometry->bufD_PeriodicRecv; @@ -1625,12 +1614,10 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, case PERIODIC_RESIDUAL: - /*--- Access the residual from the donor. ---*/ + /*--- Add contributions to total residual. ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - Residual[iVar] = bufDRecv[buf_offset]; - buf_offset++; - } + LinSysRes.AddBlock(iPoint, &bufDRecv[buf_offset]); + buf_offset += nVar; /*--- Check the computed time step against the donor value and keep the minimum in order to be conservative. ---*/ @@ -1640,21 +1627,6 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, base_nodes->SetDelta_Time(iPoint,bufDRecv[buf_offset]); buf_offset++; - /*--- Access the Jacobian from the donor if implicit. ---*/ - - if (implicit_periodic) { - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_i[iVar][jVar] = bufDRecv[buf_offset]; - buf_offset++; - } - } - } - - /*--- Add contributions to total residual. ---*/ - - LinSysRes.AddBlock(iPoint, Residual); - /*--- For implicit integration, we choose the first periodic face of each pair to be the master/owner of the solution for the linear system while fixing the @@ -1665,6 +1637,13 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, if (implicit_periodic) { + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_i[iVar][jVar] = bufDRecv[buf_offset]; + buf_offset++; + } + } + Jacobian.AddBlock2Diag(iPoint, Jacobian_i); if (iPeriodic == val_periodic_index + nPeriodic/2) { @@ -1691,19 +1670,13 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, if ((implicit_periodic) && (iPeriodic == val_periodic_index + nPeriodic/2)) { - /*--- Access the solution from the donor. ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Solution[iVar] = bufDRecv[buf_offset]; - buf_offset++; - } - /*--- Directly set the solution on the passive periodic face that is provided from the master. ---*/ for (iVar = 0; iVar < nVar; iVar++) { - base_nodes->SetSolution(iPoint, iVar, Solution[iVar]); - base_nodes->SetSolution_Old(iPoint, iVar, Solution[iVar]); + base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset]); + base_nodes->SetSolution_Old(iPoint, iVar, bufDRecv[buf_offset]); + buf_offset++; } } @@ -1900,6 +1873,11 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, delete [] Diff; + if (Jacobian_i) + for (iVar = 0; iVar < nVar; iVar++) + delete [] Jacobian_i[iVar]; + delete [] Jacobian_i; + } void CSolver::InitiateComms(CGeometry *geometry, @@ -2727,80 +2705,57 @@ void CSolver::SetResidual_BGS(CGeometry *geometry, CConfig *config) { void CSolver::SetRotatingFrame_GCL(CGeometry *geometry, CConfig *config) { - unsigned short iDim, nDim = geometry->GetnDim(), iVar, nVar = GetnVar(), iMarker; - unsigned long iVertex, iEdge; - su2double ProjGridVel; - const su2double* Normal; + /*--- Loop interior points ---*/ - /*--- Loop interior edges ---*/ + SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain,2*omp_get_max_threads())) + for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) { - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { + const su2double* GridVel_i = geometry->nodes->GetGridVel(iPoint); + const su2double* Solution_i = base_nodes->GetSolution(iPoint); - const unsigned long iPoint = geometry->edges->GetNode(iEdge,0); - const unsigned long jPoint = geometry->edges->GetNode(iEdge,1); + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) { - /*--- Solution at each edge point ---*/ + const auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); + const su2double* Normal = geometry->edges->GetNormal(iEdge); - su2double *Solution_i = base_nodes->GetSolution(iPoint); - su2double *Solution_j = base_nodes->GetSolution(jPoint); + const auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + const su2double* GridVel_j = geometry->nodes->GetGridVel(jPoint); - for (iVar = 0; iVar < nVar; iVar++) - Solution[iVar] = 0.5* (Solution_i[iVar] + Solution_j[iVar]); - - /*--- Grid Velocity at each edge point ---*/ - - su2double *GridVel_i = geometry->nodes->GetGridVel(iPoint); - su2double *GridVel_j = geometry->nodes->GetGridVel(jPoint); - for (iDim = 0; iDim < nDim; iDim++) - Vector[iDim] = 0.5* (GridVel_i[iDim] + GridVel_j[iDim]); + /*--- Determine whether to consider the normal outward or inward. ---*/ + su2double dir = (geometry->edges->GetNode(iEdge,0) == iPoint)? 0.5 : -0.5; - Normal = geometry->edges->GetNormal(iEdge); - - ProjGridVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjGridVel += Vector[iDim]*Normal[iDim]; - - for (iVar = 0; iVar < nVar; iVar++) - Residual[iVar] = ProjGridVel*Solution_i[iVar]; - - LinSysRes.AddBlock(iPoint, Residual); - - for (iVar = 0; iVar < nVar; iVar++) - Residual[iVar] = ProjGridVel*Solution_j[iVar]; - - LinSysRes.SubtractBlock(jPoint, Residual); + su2double Flux = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) + Flux += dir*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + for (auto iVar = 0u; iVar < nVar; iVar++) + LinSysRes(iPoint,iVar) += Flux * Solution_i[iVar]; + } } /*--- Loop boundary edges ---*/ - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { + for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); iMarker++) { if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - const unsigned long Point = geometry->vertex[iMarker][iVertex]->GetNode(); - /*--- Solution at each edge point ---*/ + SU2_OMP_FOR_STAT(32) + for (auto iVertex = 0u; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - su2double *Solution = base_nodes->GetSolution(Point); + const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); /*--- Grid Velocity at each edge point ---*/ - su2double *GridVel = geometry->nodes->GetGridVel(Point); + const su2double* GridVel = geometry->nodes->GetGridVel(iPoint); /*--- Summed normal components ---*/ - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - - ProjGridVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjGridVel += GridVel[iDim]*Normal[iDim]; - - for (iVar = 0; iVar < nVar; iVar++) - Residual[iVar] = ProjGridVel*Solution[iVar]; + const su2double* Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - LinSysRes.SubtractBlock(Point, Residual); + su2double Flux = GeometryToolbox::DotProduct(nDim, Normal, GridVel); + for (auto iVar = 0u; iVar < nVar; iVar++) + LinSysRes(iPoint,iVar) -= Flux * base_nodes->GetSolution(iPoint,iVar); } } } @@ -2876,6 +2831,7 @@ void CSolver::Update_Cross_Term(CConfig *config, su2passivematrix &cross_term) { * When "alpha" is 1, i.e. no relaxation, we effectively subtract the old * value and add the new one to the total ("External"). ---*/ + vector solution(nVar); passivedouble alpha = SU2_TYPE::GetValue(config->GetAitkenStatRelax()); for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { @@ -2885,10 +2841,10 @@ void CSolver::Update_Cross_Term(CConfig *config, su2passivematrix &cross_term) { delta = alpha * (new_val - cross_term(iPoint,iVar)); /*--- Update cross term. ---*/ cross_term(iPoint,iVar) += delta; - Solution[iVar] = delta; + solution[iVar] = delta; } /*--- Update the sum of all cross-terms. ---*/ - base_nodes->Add_External(iPoint, Solution); + base_nodes->Add_External(iPoint, solution.data()); } } diff --git a/SU2_CFD/src/solvers/CTransLMSolver.cpp b/SU2_CFD/src/solvers/CTransLMSolver.cpp index 0c84248b8d97..6c0017b29b24 100644 --- a/SU2_CFD/src/solvers/CTransLMSolver.cpp +++ b/SU2_CFD/src/solvers/CTransLMSolver.cpp @@ -105,18 +105,6 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh } - /*--- Computation of gradients by least squares ---*/ - if (config->GetLeastSquaresRequired()) { - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - /*--- c vector := transpose(WA)*(Wb) ---*/ - Cvector = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - /*--- Read farfield conditions from config ---*/ Intermittency_Inf = config->GetIntermittency_FreeStream(); tu_Inf = config->GetTurbulenceIntensity_FreeStream(); diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 8aaa43915661..cc6f3bfaef58 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -28,6 +28,7 @@ #include "../../include/solvers/CTurbSASolver.hpp" #include "../../include/variables/CTurbSAVariable.hpp" #include "../../../Common/include/omp_structure.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" CTurbSASolver::CTurbSASolver(void) : CTurbSolver() { } @@ -35,7 +36,7 @@ CTurbSASolver::CTurbSASolver(void) : CTurbSolver() { } CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned short iMesh, CFluidModel* FluidModel) : CTurbSolver(geometry, config) { - unsigned short iVar, iDim, nLineLets; + unsigned short iVar, nLineLets; unsigned long iPoint; su2double Density_Inf, Viscosity_Inf, Factor_nu_Inf, Factor_nu_Engine, Factor_nu_ActDisk; @@ -65,10 +66,7 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor /*--- Define some auxiliar vector related with the residual ---*/ - Residual = new su2double[nVar](); Residual_RMS = new su2double[nVar](); - Residual_i = new su2double[nVar](); - Residual_j = new su2double[nVar](); Residual_Max = new su2double[nVar](); /*--- Define some structures for locating max residuals ---*/ @@ -79,20 +77,6 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor Point_Max_Coord[iVar] = new su2double[nDim](); } - /*--- Define some auxiliar vector related with the solution ---*/ - - Solution = new su2double[nVar]; - Solution_i = new su2double[nVar]; Solution_j = new su2double[nVar]; - - /*--- Jacobians and vector structures for implicit computations ---*/ - - Jacobian_i = new su2double* [nVar]; - Jacobian_j = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double [nVar]; - Jacobian_j[iVar] = new su2double [nVar]; - } - /*--- Initialization of the structure of the whole Jacobian ---*/ if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (SA model)." << endl; @@ -116,20 +100,6 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor OutputHeadingNames = new string[nOutputVariables]; } - /*--- Computation of gradients by least squares ---*/ - - if (config->GetLeastSquaresRequired()) { - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - - /*--- c vector := transpose(WA)*(Wb) ---*/ - Cvector = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - /*--- Initialize the BGS residuals in multizone problems. ---*/ if (multizone){ Residual_BGS = new su2double[nVar](); @@ -463,63 +433,32 @@ void CTurbSASolver::Source_Template(CGeometry *geometry, CSolver **solver_contai void CTurbSASolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned long iPoint, iVertex; - unsigned short iVar; - - /*--- The dirichlet condition is used only without wall function, otherwise the - convergence is compromised as we are providing nu tilde values for the - first point of the wall ---*/ - - if (!config->GetWall_Functions()) { - - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Get the velocity vector ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - Solution[iVar] = 0.0; - - nodes->SetSolution_Old(iPoint,Solution); - LinSysRes.SetBlock_Zero(iPoint); - /*--- Includes 1 in the diagonal ---*/ - - Jacobian.DeleteValsRowi(iPoint); - } - } - } - else { - - /*--- Evaluate nu tilde at the closest point to the surface using the wall functions ---*/ + /*--- Evaluate nu tilde at the closest point to the surface using the wall functions. ---*/ + if (config->GetWall_Functions()) { + SU2_OMP_MASTER SetNuTilde_WF(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); - + SU2_OMP_BARRIER + return; } -} + /*--- The dirichlet condition is used only without wall function, otherwise the + convergence is compromised as we are providing nu tilde values for the + first point of the wall. ---*/ -void CTurbSASolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned long iPoint, iVertex; - unsigned short iVar; + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ if (geometry->nodes->GetDomain(iPoint)) { - /*--- Get the velocity vector ---*/ - for (iVar = 0; iVar < nVar; iVar++) - Solution[iVar] = 0.0; + for (auto iVar = 0u; iVar < nVar; iVar++) + nodes->SetSolution_Old(iPoint,iVar,0.0); - nodes->SetSolution_Old(iPoint,Solution); LinSysRes.SetBlock_Zero(iPoint); /*--- Includes 1 in the diagonal ---*/ @@ -530,18 +469,20 @@ void CTurbSASolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_con } -void CTurbSASolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { +void CTurbSASolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned long iPoint, iVertex; - unsigned short iVar, iDim; - su2double *Normal, *V_infty, *V_domain; + BC_HeatFlux_Wall(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); - Normal = new su2double[nDim]; +} + +void CTurbSASolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ @@ -549,11 +490,11 @@ void CTurbSASolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container /*--- Allocate the value at the infinity ---*/ - V_infty = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + auto V_infty = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Grid Movement ---*/ @@ -564,16 +505,13 @@ void CTurbSASolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container /*--- Set turbulent variable at the wall, and at infinity ---*/ - for (iVar = 0; iVar < nVar; iVar++) - Solution_i[iVar] = nodes->GetSolution(iPoint,iVar); - Solution_j[0] = nu_tilde_Inf; - conv_numerics->SetTurbVar(Solution_i, Solution_j); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), &nu_tilde_Inf); /*--- Set Normal (it is necessary to change the sign) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) - Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); conv_numerics->SetNormal(Normal); /*--- Compute residuals and Jacobians ---*/ @@ -588,26 +526,17 @@ void CTurbSASolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container } } - delete [] Normal; - } void CTurbSASolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned short iDim; - unsigned long iVertex, iPoint; - su2double *V_inlet, *V_domain, *Normal; - - Normal = new su2double[nDim]; - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ @@ -615,30 +544,26 @@ void CTurbSASolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CN /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); /*--- Allocate the value at the inlet ---*/ - V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_inlet); /*--- Set the turbulent variable states (prescribed for an inflow) ---*/ - - Solution_i[0] = nodes->GetSolution(iPoint,0); - /*--- Load the inlet turbulence variable (uniform by default). ---*/ - Solution_j[0] = Inlet_TurbVars[val_marker][iVertex][0]; - - conv_numerics->SetTurbVar(Solution_i, Solution_j); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), Inlet_TurbVars[val_marker][iVertex]); /*--- Set various other quantities in the conv_numerics class ---*/ @@ -683,24 +608,17 @@ void CTurbSASolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CN } } - /*--- Free locally allocated memory ---*/ - delete[] Normal; - } void CTurbSASolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned long iPoint, iVertex; - unsigned short iVar, iDim; - su2double *V_outlet, *V_domain, *Normal; - - Normal = new su2double[nDim]; - /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ @@ -708,11 +626,11 @@ void CTurbSASolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, C /*--- Allocate the value at the outlet ---*/ - V_outlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + auto V_outlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ @@ -720,21 +638,15 @@ void CTurbSASolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, C /*--- Set the turbulent variables. Here we use a Neumann BC such that the turbulent variable is copied from the interior of the - domain to the outlet before computing the residual. - Solution_i --> TurbVar_internal, - Solution_j --> TurbVar_outlet ---*/ + domain to the outlet before computing the residual. ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - Solution_i[iVar] = nodes->GetSolution(iPoint,iVar); - Solution_j[iVar] = nodes->GetSolution(iPoint,iVar); - } - conv_numerics->SetTurbVar(Solution_i, Solution_j); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint)); /*--- Set Normal (negate for outward convention) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) - Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); conv_numerics->SetNormal(Normal); if (dynamic_grid) @@ -776,25 +688,17 @@ void CTurbSASolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, C } } - /*--- Free locally allocated memory ---*/ - - delete[] Normal; - } void CTurbSASolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned long iPoint, iVertex; - unsigned short iDim; - su2double *V_inflow, *V_domain, *Normal; - - Normal = new su2double[nDim]; /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ @@ -802,11 +706,11 @@ void CTurbSASolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_conta /*--- Allocate the value at the infinity ---*/ - V_inflow = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + auto V_inflow = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ @@ -820,9 +724,9 @@ void CTurbSASolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_conta /*--- Set Normal (negate for outward convention) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) - Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); conv_numerics->SetNormal(Normal); /*--- Set grid movement ---*/ @@ -867,28 +771,17 @@ void CTurbSASolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_conta } - /*--- Free locally allocated memory ---*/ - - delete[] Normal; - } void CTurbSASolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned short iDim; - unsigned long iVertex, iPoint; - su2double *V_exhaust, *V_domain, *Normal; - - Normal = new su2double[nDim]; - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ @@ -896,16 +789,17 @@ void CTurbSASolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_cont /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); /*--- Allocate the value at the infinity ---*/ - V_exhaust = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + auto V_exhaust = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ @@ -913,10 +807,7 @@ void CTurbSASolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_cont /*--- Set the turbulent variable states (prescribed for an inflow) ---*/ - Solution_i[0] = nodes->GetSolution(iPoint,0); - Solution_j[0] = nu_tilde_Engine; - - conv_numerics->SetTurbVar(Solution_i, Solution_j); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), &nu_tilde_Engine); /*--- Set various other quantities in the conv_numerics class ---*/ @@ -963,10 +854,6 @@ void CTurbSASolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_cont } } - /*--- Free locally allocated memory ---*/ - - delete[] Normal; - } void CTurbSASolver::BC_ActDisk_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, @@ -985,116 +872,104 @@ void CTurbSASolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker, bool val_inlet_surface) { - unsigned long iPoint, iVertex, GlobalIndex_donor, GlobalIndex; - su2double *V_outlet, *V_inlet, *V_domain, *Normal, *UnitNormal, Area, Vn; - bool ReverseFlow; - unsigned short iDim; - - Normal = new su2double[nDim]; - UnitNormal = new su2double[nDim]; - /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - GlobalIndex_donor = solver_container[FLOW_SOL]->GetDonorGlobalIndex(val_marker, iVertex); - GlobalIndex = geometry->nodes->GetGlobalIndex(iPoint); + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto GlobalIndex_donor = solver_container[FLOW_SOL]->GetDonorGlobalIndex(val_marker, iVertex); + const auto GlobalIndex = geometry->nodes->GetGlobalIndex(iPoint); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ - if ((geometry->nodes->GetDomain(iPoint)) && (GlobalIndex != GlobalIndex_donor)) { - - /*--- Normal vector for this vertex (negate for outward convention) ---*/ - - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; - conv_numerics->SetNormal(Normal); - - Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; - Area = sqrt (Area); + if (!geometry->nodes->GetDomain(iPoint) || (GlobalIndex == GlobalIndex_donor)) { + continue; + } - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = Normal[iDim]/Area; + /*--- Normal vector for this vertex (negate for outward convention) ---*/ - /*--- Retrieve solution at the farfield boundary node ---*/ + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + su2double Area = GeometryToolbox::Norm(nDim, Normal); - /*--- Check the flow direction. Project the flow into the normal to the inlet face ---*/ + su2double UnitNormal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + UnitNormal[iDim] = Normal[iDim]/Area; - Vn = 0.0; ReverseFlow = false; - for (iDim = 0; iDim < nDim; iDim++) { Vn += V_domain[iDim+1]*UnitNormal[iDim]; } + /*--- Retrieve solution at the farfield boundary node ---*/ - if ((val_inlet_surface) && (Vn < 0.0)) { ReverseFlow = true; } - if ((!val_inlet_surface) && (Vn > 0.0)) { ReverseFlow = true; } + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); - /*--- Do not anything if there is a - reverse flow, Euler b.c. for the direct problem ---*/ + /*--- Check the flow direction. Project the flow into the normal to the inlet face ---*/ - if (!ReverseFlow) { + su2double Vn = GeometryToolbox::DotProduct(nDim, &V_domain[1], UnitNormal); - /*--- Allocate the value at the infinity ---*/ + bool ReverseFlow = false; + if ((val_inlet_surface) && (Vn < 0.0)) { ReverseFlow = true; } + if ((!val_inlet_surface) && (Vn > 0.0)) { ReverseFlow = true; } - if (val_inlet_surface) { - V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); - V_outlet = solver_container[FLOW_SOL]->GetDonorPrimVar(val_marker, iVertex); - conv_numerics->SetPrimitive(V_domain, V_inlet); - } - else { - V_outlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); - V_inlet = solver_container[FLOW_SOL]->GetDonorPrimVar(val_marker, iVertex); - conv_numerics->SetPrimitive(V_domain, V_outlet); - } + /*--- Do not anything if there is a + reverse flow, Euler b.c. for the direct problem ---*/ - /*--- Set the turb. variable solution - set the turbulent variables. Here we use a Neumann BC such - that the turbulent variable is copied from the interior of the - domain to the outlet before computing the residual. - or set the turbulent variable states (prescribed for an inflow) ----*/ + if (ReverseFlow) continue; - Solution_i[0] = nodes->GetSolution(iPoint,0); + /*--- Allocate the value at the infinity ---*/ - // if (val_inlet_surface) Solution_j[0] = 0.5*(nodes->GetSolution(iPoint,0)+V_outlet [nDim+9]); - // else Solution_j[0] = 0.5*(nodes->GetSolution(iPoint,0)+V_inlet [nDim+9]); + if (val_inlet_surface) { + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + conv_numerics->SetPrimitive(V_domain, V_inlet); + } + else { + auto V_outlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + conv_numerics->SetPrimitive(V_domain, V_outlet); + } - // /*--- Inflow analysis (interior extrapolation) ---*/ - // if (((val_inlet_surface) && (!ReverseFlow)) || ((!val_inlet_surface) && (ReverseFlow))) { - // Solution_j[0] = 2.0*node[iPoint]->GetSolution(0) - node[iPoint_Normal]->GetSolution(0); - // } + /*--- Set the turb. variable solution + set the turbulent variables. Here we use a Neumann BC such + that the turbulent variable is copied from the interior of the + domain to the outlet before computing the residual. + or set the turbulent variable states (prescribed for an inflow) ----*/ - // /*--- Outflow analysis ---*/ - // else { - // if (val_inlet_surface) Solution_j[0] = Factor_nu_ActDisk*V_outlet [nDim+9]; - // else { Solution_j[0] = Factor_nu_ActDisk*V_inlet [nDim+9]; } - // } + // if (val_inlet_surface) Solution_j[0] = 0.5*(nodes->GetSolution(iPoint,0)+V_outlet [nDim+9]); + // else Solution_j[0] = 0.5*(nodes->GetSolution(iPoint,0)+V_inlet [nDim+9]); - /*--- Inflow analysis (interior extrapolation) ---*/ - if (((val_inlet_surface) && (!ReverseFlow)) || ((!val_inlet_surface) && (ReverseFlow))) { - Solution_j[0] = nodes->GetSolution(iPoint,0); - } + // /*--- Inflow analysis (interior extrapolation) ---*/ + // if (((val_inlet_surface) && (!ReverseFlow)) || ((!val_inlet_surface) && (ReverseFlow))) { + // Solution_j[0] = 2.0*node[iPoint]->GetSolution(0) - node[iPoint_Normal]->GetSolution(0); + // } - /*--- Outflow analysis ---*/ - else { - Solution_j[0] = nu_tilde_ActDisk; - } + // /*--- Outflow analysis ---*/ + // else { + // if (val_inlet_surface) Solution_j[0] = Factor_nu_ActDisk*V_outlet [nDim+9]; + // else { Solution_j[0] = Factor_nu_ActDisk*V_inlet [nDim+9]; } + // } - conv_numerics->SetTurbVar(Solution_i, Solution_j); + if (((val_inlet_surface) && (!ReverseFlow)) || ((!val_inlet_surface) && (ReverseFlow))) { + /*--- Inflow analysis (interior extrapolation) ---*/ + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint)); + } + else { + /*--- Outflow analysis ---*/ + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), &nu_tilde_ActDisk); + } - /*--- Grid Movement ---*/ + /*--- Grid Movement ---*/ - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); - /*--- Compute the residual using an upwind scheme ---*/ + /*--- Compute the residual using an upwind scheme ---*/ - auto residual = conv_numerics->ComputeResidual(config); - LinSysRes.AddBlock(iPoint, residual); + auto residual = conv_numerics->ComputeResidual(config); + LinSysRes.AddBlock(iPoint, residual); - /*--- Jacobian contribution for implicit integration ---*/ + /*--- Jacobian contribution for implicit integration ---*/ - Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); // /*--- Viscous contribution, commented out because serious convergence problems ---*/ // @@ -1121,57 +996,47 @@ void CTurbSASolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, // LinSysRes.SubtractBlock(iPoint, residual); // Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); - } - } } - /*--- Free locally allocated memory ---*/ - - delete[] Normal; - delete[] UnitNormal; - } void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned short iDim, iSpan; - unsigned long oldVertex, iPoint, Point_Normal, iVertex; - su2double *V_inlet, *V_domain, *Normal; - su2double extAverageNu; - Normal = new su2double[nDim]; - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - unsigned short nSpanWiseSections = config->GetnSpanWiseSections(); + const auto nSpanWiseSections = config->GetnSpanWiseSections(); /*--- Loop over all the vertices on this boundary marker ---*/ - for (iSpan= 0; iSpan < nSpanWiseSections ; iSpan++){ - extAverageNu = solver_container[FLOW_SOL]->GetExtAverageNu(val_marker, iSpan); + for (auto iSpan = 0u; iSpan < nSpanWiseSections ; iSpan++){ + + su2double extAverageNu = solver_container[FLOW_SOL]->GetExtAverageNu(val_marker, iSpan); /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { /*--- find the node related to the vertex ---*/ - iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); + const auto iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); /*--- using the other vertex information for retrieving some information ---*/ - oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); + const auto oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); /*--- Index of the closest interior node ---*/ - Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); + const auto Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][oldVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][oldVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); /*--- Allocate the value at the inlet ---*/ - V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, oldVertex); + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, oldVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ @@ -1179,24 +1044,13 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c /*--- Set the turbulent variable states (prescribed for an inflow) ---*/ - Solution_i[0] = nodes->GetSolution(iPoint,0); - Solution_j[0] = extAverageNu; - - conv_numerics->SetTurbVar(Solution_i, Solution_j); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), &extAverageNu); /*--- Set various other quantities in the conv_numerics class ---*/ - conv_numerics->SetNormal(Normal); - - conv_numerics->SetTurbVar(Solution_i, Solution_j); - - /*--- Set various other quantities in the conv_numerics class ---*/ - - conv_numerics->SetNormal(Normal); - if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), - geometry->nodes->GetGridVel(iPoint)); + geometry->nodes->GetGridVel(iPoint)); /*--- Compute the residual using an upwind scheme ---*/ @@ -1209,7 +1063,8 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c /*--- Viscous contribution ---*/ - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(Point_Normal)); + visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), + geometry->nodes->GetCoord(Point_Normal)); visc_numerics->SetNormal(Normal); /*--- Conservative variables w/o reconstruction ---*/ @@ -1218,8 +1073,10 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c /*--- Turbulent variables w/o reconstruction, and its gradients ---*/ - visc_numerics->SetTurbVar(Solution_i, Solution_j); - visc_numerics->SetTurbVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(iPoint)); + visc_numerics->SetTurbVar(nodes->GetSolution(iPoint), &extAverageNu); + + visc_numerics->SetTurbVarGradient(nodes->GetGradient(iPoint), + nodes->GetGradient(iPoint)); /*--- Compute residual, and Jacobians ---*/ @@ -1233,61 +1090,55 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c } } - /*--- Free locally allocated memory ---*/ - delete[] Normal; - } void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned short iDim, iSpan; - unsigned long oldVertex, iPoint, Point_Normal, iVertex; - su2double *V_inlet, *V_domain, *Normal; - - su2double rho, pressure, muLam, Factor_nu_Inf, nu_tilde; - Normal = new su2double[nDim]; - unsigned short nSpanWiseSections = config->GetnSpanWiseSections(); - CFluidModel *FluidModel; + const auto nSpanWiseSections = config->GetnSpanWiseSections(); - FluidModel = solver_container[FLOW_SOL]->GetFluidModel(); - Factor_nu_Inf = config->GetNuFactor_FreeStream(); + CFluidModel *FluidModel = solver_container[FLOW_SOL]->GetFluidModel(); + su2double Factor_nu_Inf = config->GetNuFactor_FreeStream(); /*--- Loop over all the spans on this boundary marker ---*/ - for (iSpan= 0; iSpan < nSpanWiseSections ; iSpan++){ - rho = solver_container[FLOW_SOL]->GetAverageDensity(val_marker, iSpan); - pressure = solver_container[FLOW_SOL]->GetAveragePressure(val_marker, iSpan); + for (auto iSpan = 0; iSpan < nSpanWiseSections; iSpan++) { - FluidModel->SetTDState_Prho(pressure, rho); - muLam = FluidModel->GetLaminarViscosity(); + su2double rho = solver_container[FLOW_SOL]->GetAverageDensity(val_marker, iSpan); + su2double pressure = solver_container[FLOW_SOL]->GetAveragePressure(val_marker, iSpan); - nu_tilde = Factor_nu_Inf*muLam/rho; + FluidModel->SetTDState_Prho(pressure, rho); + su2double muLam = FluidModel->GetLaminarViscosity(); + su2double nu_tilde = Factor_nu_Inf*muLam/rho; /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { /*--- find the node related to the vertex ---*/ - iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); + const auto iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); /*--- using the other vertex information for retrieving some information ---*/ - oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); + const auto oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); /*--- Index of the closest interior node ---*/ - Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); + const auto Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][oldVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][oldVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); /*--- Allocate the value at the inlet ---*/ - V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, oldVertex); + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, oldVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ @@ -1295,24 +1146,11 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain /*--- Set the turbulent variable states (prescribed for an inflow) ---*/ - Solution_i[0] = nodes->GetSolution(iPoint,0); - Solution_j[0] = nu_tilde; - - conv_numerics->SetTurbVar(Solution_i, Solution_j); - - /*--- Set various other quantities in the conv_numerics class ---*/ - - conv_numerics->SetNormal(Normal); - - conv_numerics->SetTurbVar(Solution_i, Solution_j); - - /*--- Set various other quantities in the conv_numerics class ---*/ - - conv_numerics->SetNormal(Normal); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), &nu_tilde); if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), - geometry->nodes->GetGridVel(iPoint)); + geometry->nodes->GetGridVel(iPoint)); /*--- Compute the residual using an upwind scheme ---*/ @@ -1325,7 +1163,9 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain /*--- Viscous contribution ---*/ - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(Point_Normal)); + visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), + geometry->nodes->GetCoord(Point_Normal)); + visc_numerics->SetNormal(Normal); /*--- Conservative variables w/o reconstruction ---*/ @@ -1334,8 +1174,10 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain /*--- Turbulent variables w/o reconstruction, and its gradients ---*/ - visc_numerics->SetTurbVar(Solution_i, Solution_j); - visc_numerics->SetTurbVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(iPoint)); + visc_numerics->SetTurbVar(nodes->GetSolution(iPoint), &nu_tilde); + + visc_numerics->SetTurbVarGradient(nodes->GetGradient(iPoint), + nodes->GetGradient(iPoint)); /*--- Compute residual, and Jacobians ---*/ @@ -1349,9 +1191,6 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain } } - /*--- Free locally allocated memory ---*/ - delete[] Normal; - } void CTurbSASolver::BC_Interface_Boundary(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, @@ -1683,275 +1522,148 @@ void CTurbSASolver::BC_NearField_Boundary(CGeometry *geometry, CSolver **solver_ void CTurbSASolver::SetNuTilde_WF(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - /*--- Local variables ---*/ + const su2double Gas_Constant = config->GetGas_ConstantND(); + const su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; - unsigned short iDim, jDim, iVar, iNode; - unsigned long iVertex, iPoint, iPoint_Neighbor, counter; - - su2double func, func_prim; - su2double *Normal, Area; - su2double div_vel, UnitNormal[3]; - su2double **grad_primvar, tau[3][3]; - su2double Vel[3], VelNormal, VelTang[3], VelTangMod, VelInfMod, WallDist[3], WallDistMod; - su2double Lam_Visc_Normal, Kin_Visc_Normal, dypw_dyp, Eddy_Visc, nu_til_old, nu_til, cv1_3; - su2double T_Normal, P_Normal, Density_Normal; - su2double Density_Wall, T_Wall, P_Wall, Lam_Visc_Wall, Tau_Wall, Tau_Wall_Old; - su2double *Coord, *Coord_Normal; - su2double diff, Delta; - su2double U_Tau, U_Plus = 0.0, Gam = 0.0, Beta = 0.0, Phi, Q = 0.0, Y_Plus_White = 0.0, Y_Plus; - su2double TauElem[3], TauNormal, TauTangent[3], WallShearStress; - su2double Gas_Constant = config->GetGas_ConstantND(); - su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; - - unsigned short max_iter = 100; - su2double tol = 1e-10; - - /*--- Get the freestream velocity magnitude for non-dim. purposes ---*/ - - su2double *VelInf = config->GetVelocity_FreeStreamND(); - VelInfMod = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - VelInfMod += VelInf[iDim]; - VelInfMod = sqrt(VelInfMod); + constexpr unsigned short max_iter = 100; + const su2double tol = 1e-10; /*--- Compute the recovery factor ---*/ // su2double-check: laminar or turbulent Pr for this? - su2double Recovery = pow(config->GetPrandtl_Lam(),(1.0/3.0)); + const su2double Recovery = pow(config->GetPrandtl_Lam(),(1.0/3.0)); /*--- Typical constants from boundary layer theory ---*/ - su2double kappa = 0.4; - su2double B = 5.5; - - /*--- Identify the boundary by string name ---*/ - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - - /*--- Get the specified wall heat flux from config ---*/ + const su2double kappa = 0.4; + const su2double B = 5.5; + const su2double cv1_3 = 7.1*7.1*7.1; - // Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag); + CVariable* flow_nodes = solver_container[FLOW_SOL]->GetNodes(); /*--- Loop over all of the vertices on this boundary marker ---*/ - for(iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- We can use also GetNormal_Neighbor, and eliminate the following loop ---*/ - - iPoint_Neighbor = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - - for(iNode = 0; iNode < geometry->nodes->GetnPoint(iPoint); iNode++) { - iPoint_Neighbor = geometry->nodes->GetPoint(iPoint, iNode); - - /*--- Check if the node belongs to the domain (i.e, not a halo node) - and the neighbor is not part of the physical boundary ---*/ - - if (geometry->nodes->GetDomain(iPoint) && (!geometry->nodes->GetBoundary(iPoint_Neighbor))) { + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - /*--- Get coordinates of the current vertex and nearest normal point ---*/ + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - Coord = geometry->nodes->GetCoord(iPoint); - Coord_Normal = geometry->nodes->GetCoord(iPoint_Neighbor); + if (!geometry->nodes->GetDomain(iPoint)) continue; - /*--- Compute dual-grid area and boundary normal ---*/ + const auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); + /*--- Compute dual-grid area and boundary normal ---*/ - Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Area += Normal[iDim]*Normal[iDim]; - Area = sqrt (Area); + const auto Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = -Normal[iDim]/Area; + su2double Area = GeometryToolbox::Norm(nDim, Normal); - /*--- Get the velocity, pressure, and temperature at the nearest - (normal) interior point. ---*/ + su2double UnitNormal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + UnitNormal[iDim] = -Normal[iDim]/Area; - for (iDim = 0; iDim < nDim; iDim++) - Vel[iDim] = solver_container[FLOW_SOL]->GetNodes()->GetVelocity(iPoint_Neighbor,iDim); - P_Normal = solver_container[FLOW_SOL]->GetNodes()->GetPressure(iPoint_Neighbor); - T_Normal = solver_container[FLOW_SOL]->GetNodes()->GetTemperature(iPoint_Neighbor); + /*--- Get the velocity, pressure, and temperature at the nearest + (normal) interior point. ---*/ - /*--- Compute the wall-parallel velocity at first point off the wall ---*/ + su2double Vel[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Vel[iDim] = flow_nodes->GetVelocity(Point_Normal,iDim); + su2double P_Normal = flow_nodes->GetPressure(Point_Normal); + su2double T_Normal = flow_nodes->GetTemperature(Point_Normal); - VelNormal = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - VelNormal += Vel[iDim] * UnitNormal[iDim]; - for (iDim = 0; iDim < nDim; iDim++) - VelTang[iDim] = Vel[iDim] - VelNormal*UnitNormal[iDim]; + /*--- Compute the wall-parallel velocity at first point off the wall ---*/ - VelTangMod = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - VelTangMod += VelTang[iDim]*VelTang[iDim]; - VelTangMod = sqrt(VelTangMod); + su2double VelNormal = GeometryToolbox::DotProduct(int(MAXNDIM), Vel, UnitNormal); - /*--- Compute normal distance of the interior point from the wall ---*/ + su2double VelTang[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + VelTang[iDim] = Vel[iDim] - VelNormal*UnitNormal[iDim]; - for (iDim = 0; iDim < nDim; iDim++) - WallDist[iDim] = (Coord[iDim] - Coord_Normal[iDim]); + su2double VelTangMod = GeometryToolbox::Norm(int(MAXNDIM), VelTang); - WallDistMod = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - WallDistMod += WallDist[iDim]*WallDist[iDim]; - WallDistMod = sqrt(WallDistMod); + /*--- Compute the wall temperature using the Crocco-Buseman equation ---*/ - /*--- Compute mach number ---*/ + //T_Wall = T_Normal * (1.0 + 0.5*Gamma_Minus_One*Recovery*M_Normal*M_Normal); + su2double T_Wall = T_Normal + Recovery*pow(VelTangMod,2.0)/(2.0*Cp); - // M_Normal = VelTangMod / sqrt(Gamma * Gas_Constant * T_Normal); + /*--- Extrapolate the pressure from the interior & compute the + wall density using the equation of state ---*/ - /*--- Compute the wall temperature using the Crocco-Buseman equation ---*/ + su2double P_Wall = P_Normal; + su2double Density_Wall = P_Wall/(Gas_Constant*T_Wall); - //T_Wall = T_Normal * (1.0 + 0.5*Gamma_Minus_One*Recovery*M_Normal*M_Normal); - T_Wall = T_Normal + Recovery*pow(VelTangMod,2.0)/(2.0*Cp); + /*--- Get wall shear stress computed by flow solver ---*/ - /*--- Extrapolate the pressure from the interior & compute the - wall density using the equation of state ---*/ + su2double Lam_Visc_Wall = flow_nodes->GetLaminarViscosity(iPoint); + su2double Tau_Wall = flow_nodes->GetTauWall(iPoint); - P_Wall = P_Normal; - Density_Wall = P_Wall/(Gas_Constant*T_Wall); + /*--- Friction velocity and u+ ---*/ - /*--- Compute the shear stress at the wall in the regular fashion - by using the stress tensor on the surface ---*/ + su2double U_Tau = sqrt(Tau_Wall/Density_Wall); + su2double U_Plus = VelTangMod/U_Tau; - Lam_Visc_Wall = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); - grad_primvar = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); + /*--- Gamma, Beta, Q, and Phi, defined by Nichols & Nelson (2004) ---*/ - div_vel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - div_vel += grad_primvar[iDim+1][iDim]; + su2double Gam = Recovery*U_Tau*U_Tau/(2.0*Cp*T_Wall); + su2double Beta = 0.0; // For adiabatic flows only + su2double Q = sqrt(Beta*Beta + 4.0*Gam); + su2double Phi = asin(-1.0*Beta/Q); - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - Delta = 0.0; if (iDim == jDim) Delta = 1.0; - tau[iDim][jDim] = Lam_Visc_Wall*( grad_primvar[jDim+1][iDim] - + grad_primvar[iDim+1][jDim]) - - TWO3*Lam_Visc_Wall*div_vel*Delta; - } - TauElem[iDim] = 0.0; - for (jDim = 0; jDim < nDim; jDim++) - TauElem[iDim] += tau[iDim][jDim]*UnitNormal[jDim]; - } - - /*--- Compute wall shear stress as the magnitude of the wall-tangential - component of the shear stress tensor---*/ - - TauNormal = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - TauNormal += TauElem[iDim] * UnitNormal[iDim]; - - for (iDim = 0; iDim < nDim; iDim++) - TauTangent[iDim] = TauElem[iDim] - TauNormal * UnitNormal[iDim]; - - WallShearStress = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - WallShearStress += TauTangent[iDim]*TauTangent[iDim]; - WallShearStress = sqrt(WallShearStress); + /*--- Y+ defined by White & Christoph (compressibility and heat transfer) + negative value for (2.0*Gam*U_Plus - Beta)/Q ---*/ - /*--- Calculate the quantities from boundary layer theory and - iteratively solve for a new wall shear stress. Use the current wall - shear stress as a starting guess for the wall function. ---*/ + su2double Y_Plus_White = exp((kappa/sqrt(Gam))*(asin((2.0*Gam*U_Plus - Beta)/Q) - Phi))*exp(-1.0*kappa*B); - Tau_Wall_Old = WallShearStress; - counter = 0; diff = 1.0; + /*--- Now compute the Eddy viscosity at the first point off of the wall ---*/ - while (diff > tol) { + su2double Lam_Visc_Normal = flow_nodes->GetLaminarViscosity(Point_Normal); + su2double Density_Normal = flow_nodes->GetDensity(Point_Normal); + su2double Kin_Visc_Normal = Lam_Visc_Normal/Density_Normal; - /*--- Friction velocity and u+ ---*/ + su2double dypw_dyp = 2.0*Y_Plus_White*(kappa*sqrt(Gam)/Q)*sqrt(1.0 - pow(2.0*Gam*U_Plus - Beta,2.0)/(Q*Q)); + su2double Eddy_Visc = Lam_Visc_Wall*(1.0 + dypw_dyp - kappa*exp(-1.0*kappa*B)* + (1.0 + kappa*U_Plus + kappa*kappa*U_Plus*U_Plus/2.0) + - Lam_Visc_Normal/Lam_Visc_Wall); - U_Tau = sqrt(Tau_Wall_Old/Density_Wall); - U_Plus = VelTangMod/U_Tau; + /*--- Eddy viscosity should be always a positive number ---*/ - /*--- Gamma, Beta, Q, and Phi, defined by Nichols & Nelson (2004) ---*/ + Eddy_Visc = max(0.0, Eddy_Visc); - Gam = Recovery*U_Tau*U_Tau/(2.0*Cp*T_Wall); - Beta = 0.0; // For adiabatic flows only - Q = sqrt(Beta*Beta + 4.0*Gam); - Phi = asin(-1.0*Beta/Q); + /*--- Solve for the new value of nu_tilde given the eddy viscosity and using a Newton method ---*/ - /*--- Y+ defined by White & Christoph (compressibility and heat transfer) ---*/ + su2double nu_til_old = nodes->GetSolution(iPoint,0); - Y_Plus_White = exp((kappa/sqrt(Gam))*(asin((2.0*Gam*U_Plus - Beta)/Q) - Phi))*exp(-1.0*kappa*B); + unsigned short counter = 0; + su2double diff = 1.0; - /*--- Spalding's universal form for the BL velocity with the - outer velocity form of White & Christoph above. ---*/ + while (diff > tol) { - Y_Plus = U_Plus + Y_Plus_White - (exp(-1.0*kappa*B)* - (1.0 + kappa*U_Plus + kappa*kappa*U_Plus*U_Plus/2.0 + - kappa*kappa*kappa*U_Plus*U_Plus*U_Plus/6.0)); + su2double const_term = (Eddy_Visc/Density_Normal) * pow(Kin_Visc_Normal,3)*cv1_3; - /*--- Calculate an updated value for the wall shear stress - using the y+ value, the definition of y+, and the definition of - the friction velocity. ---*/ + su2double nu_til_2 = pow(nu_til_old,2); + su2double nu_til_3 = nu_til_old * nu_til_2; - Tau_Wall = (1.0/Density_Wall)*pow(Y_Plus*Lam_Visc_Wall/WallDistMod,2.0); + su2double func = nu_til_old * (nu_til_3 - (Eddy_Visc/Density_Normal)*nu_til_2) + const_term; + su2double func_prime = 4.0*nu_til_3 - 3.0*(Eddy_Visc/Density_Normal)*nu_til_2; + su2double nu_til = nu_til_old - func/func_prime; - /*--- Difference between the old and new Tau. Update old value. ---*/ - - diff = fabs(Tau_Wall-Tau_Wall_Old); - Tau_Wall_Old += 0.25*(Tau_Wall-Tau_Wall_Old); - - counter++; - if (counter > max_iter) { - cout << "WARNING: Tau_Wall evaluation has not converged in solver_direct_turbulent" << endl; - break; - } - - } + diff = fabs(nu_til-nu_til_old); + nu_til_old = nu_til; - /*--- Now compute the Eddy viscosity at the first point off of the wall ---*/ - - Lam_Visc_Normal = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint_Neighbor); - Density_Normal = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint_Neighbor); - Kin_Visc_Normal = Lam_Visc_Normal/Density_Normal; - - dypw_dyp = 2.0*Y_Plus_White*(kappa*sqrt(Gam)/Q)*sqrt(1.0 - pow(2.0*Gam*U_Plus - Beta,2.0)/(Q*Q)); - Eddy_Visc = Lam_Visc_Wall*(1.0 + dypw_dyp - kappa*exp(-1.0*kappa*B)* - (1.0 + kappa*U_Plus - + kappa*kappa*U_Plus*U_Plus/2.0) - - Lam_Visc_Normal/Lam_Visc_Wall); - - /*--- Eddy viscosity should be always a positive number ---*/ - - Eddy_Visc = max(0.0, Eddy_Visc); - - /*--- Solve for the new value of nu_tilde given the eddy viscosity and using a Newton method ---*/ - - nu_til_old = 0.0; nu_til = 0.0; cv1_3 = 7.1*7.1*7.1; - nu_til_old = nodes->GetSolution(iPoint,0); - counter = 0; diff = 1.0; - - while (diff > tol) { - - func = nu_til_old*nu_til_old*nu_til_old*nu_til_old - (Eddy_Visc/Density_Normal)*(nu_til_old*nu_til_old*nu_til_old + Kin_Visc_Normal*Kin_Visc_Normal*Kin_Visc_Normal*cv1_3); - func_prim = 4.0 * nu_til_old*nu_til_old*nu_til_old - 3.0*(Eddy_Visc/Density_Normal)*(nu_til_old*nu_til_old); - nu_til = nu_til_old - func/func_prim; - - diff = fabs(nu_til-nu_til_old); - nu_til_old = nu_til; - - counter++; - if (counter > max_iter) { - cout << "WARNING: Nu_tilde evaluation has not converged." << endl; - break; - } - - } - - for (iVar = 0; iVar < nVar; iVar++) - Solution[iVar] = nu_til; - - nodes->SetSolution_Old(iPoint_Neighbor,Solution); - LinSysRes.SetBlock_Zero(iPoint_Neighbor); + counter++; + if (counter > max_iter) { + cout << "WARNING: Nu_tilde evaluation has not converged." << endl; + break; + } + } - /*--- includes 1 in the diagonal ---*/ + nodes->SetSolution_Old(Point_Normal, 0, nu_til_old); + LinSysRes.SetBlock_Zero(Point_Normal); - Jacobian.DeleteValsRowi(iPoint_Neighbor); + /*--- includes 1 in the diagonal ---*/ - } + Jacobian.DeleteValsRowi(Point_Normal); - } } + } void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CConfig *config){ diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index b3e5591fc917..cfa5e773102b 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -28,13 +28,14 @@ #include "../../include/solvers/CTurbSSTSolver.hpp" #include "../../include/variables/CTurbSSTVariable.hpp" #include "../../../Common/include/omp_structure.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" CTurbSSTSolver::CTurbSSTSolver(void) : CTurbSolver() { } CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) : CTurbSolver(geometry, config) { - unsigned short iVar, iDim, nLineLets; + unsigned short iVar, nLineLets; unsigned long iPoint; ifstream restart_file; string text_line; @@ -67,10 +68,7 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- Define some auxiliary vector related with the residual ---*/ - Residual = new su2double[nVar](); Residual_RMS = new su2double[nVar](); - Residual_i = new su2double[nVar](); - Residual_j = new su2double[nVar](); Residual_Max = new su2double[nVar](); /*--- Define some structures for locating max residuals ---*/ @@ -81,20 +79,6 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh Point_Max_Coord[iVar] = new su2double[nDim](); } - /*--- Define some auxiliary vector related with the solution ---*/ - - Solution = new su2double[nVar]; - Solution_i = new su2double[nVar]; Solution_j = new su2double[nVar]; - - /*--- Jacobians and vector structures for implicit computations ---*/ - - Jacobian_i = new su2double* [nVar]; - Jacobian_j = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = new su2double [nVar]; - Jacobian_j[iVar] = new su2double [nVar]; - } - /*--- Initialization of the structure of the whole Jacobian ---*/ if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (SST model)." << endl; @@ -127,19 +111,6 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh } - /*--- Computation of gradients by least squares ---*/ - - if (config->GetLeastSquaresRequired()) { - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - Smatrix = new su2double* [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - /*--- c vector := transpose(WA)*(Wb) ---*/ - Cvector = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) - Cvector[iVar] = new su2double [nDim]; - } - /*--- Initialize value for model constants ---*/ constants[0] = 0.85; //sigma_k1 constants[1] = 1.0; //sigma_k2 @@ -168,12 +139,9 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh Intensity = config->GetTurbulenceIntensity_FreeStream(); viscRatio = config->GetTurb2LamViscRatio_FreeStream(); - su2double VelMag = 0; - for (iDim = 0; iDim < nDim; iDim++) - VelMag += VelInf[iDim]*VelInf[iDim]; - VelMag = sqrt(VelMag); + su2double VelMag2 = GeometryToolbox::SquaredNorm(nDim, VelInf); - kine_Inf = 3.0/2.0*(VelMag*VelMag*Intensity*Intensity); + kine_Inf = 3.0/2.0*(VelMag2*Intensity*Intensity); omega_Inf = rhoInf*kine_Inf/(muLamInf*viscRatio); /*--- Eddy viscosity, initialized without stress limiter at the infinity ---*/ @@ -423,46 +391,39 @@ void CTurbSSTSolver::Source_Template(CGeometry *geometry, CSolver **solver_conta void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned long iPoint, jPoint, iVertex, total_index; - unsigned short iDim, iVar; - su2double distance, density = 0.0, laminar_viscosity = 0.0, beta_1; + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ if (geometry->nodes->GetDomain(iPoint)) { /*--- distance to closest neighbor ---*/ - jPoint = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - distance = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - distance += pow(geometry->nodes->GetCoord(iPoint, iDim)- - geometry->nodes->GetCoord(jPoint, iDim), 2); - } - distance = sqrt(distance); + const auto jPoint = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + su2double distance2 = GeometryToolbox::SquaredDistance(nDim, + geometry->nodes->GetCoord(iPoint), + geometry->nodes->GetCoord(jPoint)); /*--- Set wall values ---*/ - density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(jPoint); - laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(jPoint); + su2double density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(jPoint); + su2double laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(jPoint); - beta_1 = constants[4]; + su2double beta_1 = constants[4]; - Solution[0] = 0.0; - Solution[1] = 60.0*laminar_viscosity/(density*beta_1*distance*distance); + su2double solution[2]; + solution[0] = 0.0; + solution[1] = 60.0*laminar_viscosity/(density*beta_1*distance2); /*--- Set the solution values and zero the residual ---*/ - nodes->SetSolution_Old(iPoint,Solution); - nodes->SetSolution(iPoint,Solution); + nodes->SetSolution_Old(iPoint,solution); + nodes->SetSolution(iPoint,solution); LinSysRes.SetBlock_Zero(iPoint); /*--- Change rows of the Jacobian (includes 1 in the diagonal) ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } - + Jacobian.DeleteValsRowi(iPoint*nVar); + Jacobian.DeleteValsRowi(iPoint*nVar+1); } } @@ -471,63 +432,17 @@ void CTurbSSTSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_cont void CTurbSSTSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned long iPoint, jPoint, iVertex, total_index; - unsigned short iDim, iVar; - su2double distance, density = 0.0, laminar_viscosity = 0.0, beta_1; - - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- distance to closest neighbor ---*/ - jPoint = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - distance = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - distance += (geometry->nodes->GetCoord(iPoint, iDim) - geometry->nodes->GetCoord(jPoint, iDim))* - (geometry->nodes->GetCoord(iPoint, iDim) - geometry->nodes->GetCoord(jPoint, iDim)); - } - distance = sqrt(distance); - - /*--- Set wall values ---*/ - - density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(jPoint); - laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(jPoint); - - beta_1 = constants[4]; - - Solution[0] = 0.0; - Solution[1] = 60.0*laminar_viscosity/(density*beta_1*distance*distance); - - /*--- Set the solution values and zero the residual ---*/ - nodes->SetSolution_Old(iPoint,Solution); - nodes->SetSolution(iPoint,Solution); - LinSysRes.SetBlock_Zero(iPoint); - - /*--- Change rows of the Jacobian (includes 1 in the diagonal) ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } - - } - } + BC_HeatFlux_Wall(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); } void CTurbSSTSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned long iPoint, iVertex; - su2double *Normal, *V_infty, *V_domain; - unsigned short iVar, iDim; - - Normal = new su2double[nDim]; + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ @@ -535,35 +450,32 @@ void CTurbSSTSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_containe /*--- Allocate the value at the infinity ---*/ - V_infty = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + auto V_infty = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); conv_numerics->SetPrimitive(V_domain, V_infty); /*--- Set turbulent variable at the wall, and at infinity ---*/ - for (iVar = 0; iVar < nVar; iVar++) - Solution_i[iVar] = nodes->GetSolution(iPoint,iVar); - - Solution_j[0] = kine_Inf; - Solution_j[1] = omega_Inf; + su2double solution_j[] = {kine_Inf, omega_Inf}; - conv_numerics->SetTurbVar(Solution_i, Solution_j); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), solution_j); /*--- Set Normal (it is necessary to change the sign) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) - Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); conv_numerics->SetNormal(Normal); /*--- Grid Movement ---*/ if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), + geometry->nodes->GetGridVel(iPoint)); /*--- Compute residuals and Jacobians ---*/ @@ -573,30 +485,20 @@ void CTurbSSTSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_containe LinSysRes.AddBlock(iPoint, residual); Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - } } - delete [] Normal; - } void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned short iVar, iDim; - unsigned long iVertex, iPoint; - su2double *V_inlet, *V_domain, *Normal; - - Normal = new su2double[nDim]; - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ @@ -604,16 +506,18 @@ void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); /*--- Allocate the value at the inlet ---*/ - V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ @@ -621,21 +525,13 @@ void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C /*--- Set the turbulent variable states. Use free-stream SST values for the turbulent state at the inflow. ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - Solution_i[iVar] = nodes->GetSolution(iPoint,iVar); - /*--- Load the inlet turbulence variables (uniform by default). ---*/ - Solution_j[0] = Inlet_TurbVars[val_marker][iVertex][0]; - Solution_j[1] = Inlet_TurbVars[val_marker][iVertex][1]; - - conv_numerics->SetTurbVar(Solution_i, Solution_j); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), + Inlet_TurbVars[val_marker][iVertex]); /*--- Set various other quantities in the solver class ---*/ - conv_numerics->SetNormal(Normal); - if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); @@ -680,25 +576,17 @@ void CTurbSSTSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, C } - /*--- Free locally allocated memory ---*/ - - delete [] Normal; - } void CTurbSSTSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned long iPoint, iVertex; - unsigned short iVar, iDim; - su2double *V_outlet, *V_domain, *Normal; - - Normal = new su2double[nDim]; - /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ @@ -706,11 +594,11 @@ void CTurbSSTSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, /*--- Allocate the value at the outlet ---*/ - V_outlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); + auto V_outlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, iVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ @@ -718,21 +606,16 @@ void CTurbSSTSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, /*--- Set the turbulent variables. Here we use a Neumann BC such that the turbulent variable is copied from the interior of the - domain to the outlet before computing the residual. - Solution_i --> TurbVar_internal, - Solution_j --> TurbVar_outlet ---*/ + domain to the outlet before computing the residual. ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - Solution_i[iVar] = nodes->GetSolution(iPoint,iVar); - Solution_j[iVar] = nodes->GetSolution(iPoint,iVar); - } - conv_numerics->SetTurbVar(Solution_i, Solution_j); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), + nodes->GetSolution(iPoint)); /*--- Set Normal (negate for outward convention) ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) - Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); conv_numerics->SetNormal(Normal); if (dynamic_grid) @@ -778,55 +661,49 @@ void CTurbSSTSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, } } - /*--- Free locally allocated memory ---*/ - delete[] Normal; - } void CTurbSSTSolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned short iVar, iSpan, iDim; - unsigned long oldVertex, iPoint, Point_Normal, iVertex; - su2double *V_inlet, *V_domain, *Normal; - su2double extAverageKine, extAverageOmega; - unsigned short nSpanWiseSections = config->GetnSpanWiseSections(); - - Normal = new su2double[nDim]; - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + const auto nSpanWiseSections = config->GetnSpanWiseSections(); /*--- Loop over all the vertices on this boundary marker ---*/ - for (iSpan= 0; iSpan < nSpanWiseSections ; iSpan++){ - extAverageKine = solver_container[FLOW_SOL]->GetExtAverageKine(val_marker, iSpan); - extAverageOmega = solver_container[FLOW_SOL]->GetExtAverageOmega(val_marker, iSpan); + for (auto iSpan = 0u; iSpan < nSpanWiseSections ; iSpan++){ + + su2double extAverageKine = solver_container[FLOW_SOL]->GetExtAverageKine(val_marker, iSpan); + su2double extAverageOmega = solver_container[FLOW_SOL]->GetExtAverageOmega(val_marker, iSpan); + su2double solution_j[] = {extAverageKine, extAverageOmega}; /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { /*--- find the node related to the vertex ---*/ - iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); + const auto iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); /*--- using the other vertex information for retrieving some information ---*/ - oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); + const auto oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); /*--- Index of the closest interior node ---*/ - Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); + const auto Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][oldVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][oldVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); /*--- Allocate the value at the inlet ---*/ - V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, oldVertex); + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, oldVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ @@ -834,20 +711,11 @@ void CTurbSSTSolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_ /*--- Set the turbulent variable states (prescribed for an inflow) ---*/ - for (iVar = 0; iVar < nVar; iVar++) - Solution_i[iVar] = nodes->GetSolution(iPoint,iVar); - - Solution_j[0]= extAverageKine; - Solution_j[1]= extAverageOmega; - - conv_numerics->SetTurbVar(Solution_i, Solution_j); - - /*--- Set various other quantities in the solver class ---*/ - conv_numerics->SetNormal(Normal); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), solution_j); if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), - geometry->nodes->GetGridVel(iPoint)); + geometry->nodes->GetGridVel(iPoint)); /*--- Compute the residual using an upwind scheme ---*/ auto conv_residual = conv_numerics->ComputeResidual(config); @@ -864,7 +732,7 @@ void CTurbSSTSolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_ visc_numerics->SetPrimitive(V_domain, V_inlet); /*--- Turbulent variables w/o reconstruction, and its gradients ---*/ - visc_numerics->SetTurbVar(Solution_i, Solution_j); + visc_numerics->SetTurbVar(nodes->GetSolution(iPoint), solution_j); visc_numerics->SetTurbVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(iPoint)); /*--- Menter's first blending function ---*/ @@ -880,95 +748,74 @@ void CTurbSSTSolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_ } } - /*--- Free locally allocated memory ---*/ - delete[] Normal; - } void CTurbSSTSolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - unsigned short iVar, iSpan, iDim; - unsigned long oldVertex, iPoint, Point_Normal, iVertex; - su2double *V_inlet, *V_domain, *Normal; - unsigned short nSpanWiseSections = config->GetnSpanWiseSections(); + const auto nSpanWiseSections = config->GetnSpanWiseSections(); /*--- Quantities for computing the kine and omega to impose at the inlet boundary. ---*/ - su2double rho, pressure, *Vel, VelMag, muLam, Intensity, viscRatio, kine_b, omega_b, kine; - CFluidModel *FluidModel; + CFluidModel *FluidModel = solver_container[FLOW_SOL]->GetFluidModel(); - FluidModel = solver_container[FLOW_SOL]->GetFluidModel(); - Intensity = config->GetTurbulenceIntensity_FreeStream(); - viscRatio = config->GetTurb2LamViscRatio_FreeStream(); - - Normal = new su2double[nDim]; - Vel = new su2double[nDim]; - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + su2double Intensity = config->GetTurbulenceIntensity_FreeStream(); + su2double viscRatio = config->GetTurb2LamViscRatio_FreeStream(); - - for (iSpan= 0; iSpan < nSpanWiseSections ; iSpan++){ + for (auto iSpan = 0u; iSpan < nSpanWiseSections ; iSpan++){ /*--- Compute the inflow kine and omega using the span wise averge quntities---*/ - for (iDim = 0; iDim < nDim; iDim++) - Vel[iDim] = solver_container[FLOW_SOL]->GetAverageTurboVelocity(val_marker, iSpan)[iDim]; - rho = solver_container[FLOW_SOL]->GetAverageDensity(val_marker, iSpan); - pressure = solver_container[FLOW_SOL]->GetAveragePressure(val_marker, iSpan); - kine = solver_container[FLOW_SOL]->GetAverageKine(val_marker, iSpan); + su2double rho = solver_container[FLOW_SOL]->GetAverageDensity(val_marker, iSpan); + su2double pressure = solver_container[FLOW_SOL]->GetAveragePressure(val_marker, iSpan); + su2double kine = solver_container[FLOW_SOL]->GetAverageKine(val_marker, iSpan); FluidModel->SetTDState_Prho(pressure, rho); - muLam = FluidModel->GetLaminarViscosity(); + su2double muLam = FluidModel->GetLaminarViscosity(); + + su2double VelMag2 = GeometryToolbox::SquaredNorm(nDim, + solver_container[FLOW_SOL]->GetAverageTurboVelocity(val_marker, iSpan)); - VelMag = 0; - for (iDim = 0; iDim < nDim; iDim++) - VelMag += Vel[iDim]*Vel[iDim]; - VelMag = sqrt(VelMag); + su2double kine_b = 3.0/2.0*(VelMag2*Intensity*Intensity); + su2double omega_b = rho*kine/(muLam*viscRatio); - kine_b = 3.0/2.0*(VelMag*VelMag*Intensity*Intensity); - omega_b = rho*kine/(muLam*viscRatio); + su2double solution_j[] = {kine_b, omega_b}; /*--- Loop over all the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0u; iVertex < geometry->GetnVertexSpan(val_marker,iSpan); iVertex++) { /*--- find the node related to the vertex ---*/ - iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); + const auto iPoint = geometry->turbovertex[val_marker][iSpan][iVertex]->GetNode(); /*--- using the other vertex information for retrieving some information ---*/ - oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); + const auto oldVertex = geometry->turbovertex[val_marker][iSpan][iVertex]->GetOldVertex(); /*--- Index of the closest interior node ---*/ - Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); + const auto Point_Normal = geometry->vertex[val_marker][oldVertex]->GetNormal_Neighbor(); /*--- Normal vector for this vertex (negate for outward convention) ---*/ - geometry->vertex[val_marker][oldVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + su2double Normal[MAXNDIM] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) + Normal[iDim] = -geometry->vertex[val_marker][oldVertex]->GetNormal(iDim); + conv_numerics->SetNormal(Normal); /*--- Allocate the value at the inlet ---*/ - V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, oldVertex); + auto V_inlet = solver_container[FLOW_SOL]->GetCharacPrimVar(val_marker, oldVertex); /*--- Retrieve solution at the farfield boundary node ---*/ - V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto V_domain = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_inlet); - for (iVar = 0; iVar < nVar; iVar++) - Solution_i[iVar] = nodes->GetSolution(iPoint,iVar); - /*--- Set the turbulent variable states. Use average span-wise values values for the turbulent state at the inflow. ---*/ - Solution_j[0]= kine_b; - Solution_j[1]= omega_b; - - conv_numerics->SetTurbVar(Solution_i, Solution_j); - - /*--- Set various other quantities in the solver class ---*/ - conv_numerics->SetNormal(Normal); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), solution_j); if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), @@ -982,14 +829,16 @@ void CTurbSSTSolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contai Jacobian.AddBlock2Diag(iPoint, conv_residual.jacobian_i); /*--- Viscous contribution ---*/ - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(Point_Normal)); + visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), + geometry->nodes->GetCoord(Point_Normal)); visc_numerics->SetNormal(Normal); /*--- Conservative variables w/o reconstruction ---*/ visc_numerics->SetPrimitive(V_domain, V_inlet); /*--- Turbulent variables w/o reconstruction, and its gradients ---*/ - visc_numerics->SetTurbVar(Solution_i, Solution_j); + visc_numerics->SetTurbVar(nodes->GetSolution(iPoint), solution_j); + visc_numerics->SetTurbVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(iPoint)); /*--- Menter's first blending function ---*/ @@ -1005,10 +854,6 @@ void CTurbSSTSolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contai } } - /*--- Free locally allocated memory ---*/ - delete[] Normal; - delete[] Vel; - } void CTurbSSTSolver::SetInletAtVertex(su2double *val_inlet, diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index 10444ae32e9b..e8084b961ec6 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -25,7 +25,6 @@ * License along with SU2. If not, see . */ - #include "../../include/solvers/CTurbSolver.hpp" #include "../../../Common/include/omp_structure.hpp" @@ -397,10 +396,12 @@ void CTurbSolver::BC_Periodic(CGeometry *geometry, CSolver **solver_container, accumulated corectly during the communications. For implicit calculations the Jacobians and linear system are also correctly adjusted here. ---*/ + SU2_OMP_MASTER for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_RESIDUAL); CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_RESIDUAL); } + SU2_OMP_BARRIER } @@ -410,11 +411,13 @@ void CTurbSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_conta const bool sst = (config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST); const auto nPrimVar = solver_container[FLOW_SOL]->GetnPrimVar(); su2double *PrimVar_j = new su2double[nPrimVar]; + su2double solution_j[MAXNVAR] = {0.0}; for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (config->GetMarker_All_KindBC(iMarker) != FLUID_INTERFACE) continue; + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) for (auto iVertex = 0u; iVertex < geometry->nVertex[iMarker]; iVertex++) { const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); @@ -428,16 +431,10 @@ void CTurbSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_conta for (auto iDim = 0u; iDim < nDim; iDim++) Normal[iDim] = -geometry->vertex[iMarker][iVertex]->GetNormal()[iDim]; - /*--- Initialize Residual, this will serve to accumulate the average ---*/ - - for (auto iVar = 0u; iVar < nVar; iVar++) { - Residual[iVar] = 0.0; - for (auto jVar = 0u; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - } - su2double* PrimVar_i = solver_container[FLOW_SOL]->GetNodes()->GetPrimitive(iPoint); + auto Jacobian_i = Jacobian.GetBlock(iPoint,iPoint); + /*--- Loop over the nDonorVertexes and compute the averaged flux ---*/ for (auto jVertex = 0; jVertex < nDonorVertex; jVertex++) { @@ -456,9 +453,9 @@ void CTurbSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_conta /*--- Set the turbulent variable states ---*/ for (auto iVar = 0u; iVar < nVar; ++iVar) - Solution_j[iVar] = GetSlidingState(iMarker, iVertex, iVar, jVertex); + solution_j[iVar] = GetSlidingState(iMarker, iVertex, iVar, jVertex); - conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), Solution_j); + conv_numerics->SetTurbVar(nodes->GetSolution(iPoint), solution_j); /*--- Set the normal vector ---*/ @@ -472,18 +469,12 @@ void CTurbSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_conta /*--- Accumulate the residuals to compute the average ---*/ for (auto iVar = 0u; iVar < nVar; iVar++) { - Residual[iVar] += weight*residual.residual[iVar]; + LinSysRes(iPoint,iVar) += weight*residual.residual[iVar]; for (auto jVar = 0u; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] += weight*residual.jacobian_i[iVar][jVar]; + Jacobian_i[iVar*nVar+jVar] += SU2_TYPE::GetValue(weight*residual.jacobian_i[iVar][jVar]); } } - /*--- Add Residuals and Jacobians ---*/ - - LinSysRes.AddBlock(iPoint, Residual); - - Jacobian.AddBlock2Diag(iPoint, Jacobian_i); - /*--- Set the normal vector and the coordinates ---*/ visc_numerics->SetNormal(Normal); @@ -495,7 +486,7 @@ void CTurbSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_conta /*--- Turbulent variables and their gradients ---*/ - visc_numerics->SetTurbVar(Solution_i, Solution_j); + visc_numerics->SetTurbVar(nodes->GetSolution(iPoint), solution_j); visc_numerics->SetTurbVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(iPoint)); /*--- Menter's first blending function ---*/ diff --git a/SU2_PY/pySU2/pySU2.i b/SU2_PY/pySU2/pySU2.i index 546c4aba6c88..31fd5d4cabb8 100644 --- a/SU2_PY/pySU2/pySU2.i +++ b/SU2_PY/pySU2/pySU2.i @@ -46,7 +46,7 @@ threads="1" %} // ----------- USED MODULES ------------ -%import "../../Common/include/datatype_structure.hpp" +%import "../../Common/include/basic_types/datatype_structure.hpp" %import "../../Common/include/mpi_structure.hpp" %include "std_string.i" %include "std_vector.i" diff --git a/SU2_PY/pySU2/pySU2ad.i b/SU2_PY/pySU2/pySU2ad.i index 76ad6f0cf58b..d2b5c92ce98f 100644 --- a/SU2_PY/pySU2/pySU2ad.i +++ b/SU2_PY/pySU2/pySU2ad.i @@ -46,7 +46,7 @@ threads="1" %} // ----------- USED MODULES ------------ -%import "../../Common/include/datatype_structure.hpp" +%import "../../Common/include/basic_types/datatype_structure.hpp" %import "../../Common/include/mpi_structure.hpp" %include "std_string.i" %include "std_vector.i" diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py new file mode 100644 index 000000000000..5fef6d079c4b --- /dev/null +++ b/TestCases/hybrid_regression.py @@ -0,0 +1,594 @@ +#!/usr/bin/env python + +## \file parallel_regression.py +# \brief Python script for automated regression testing of SU2 examples +# \author A. Aranake, A. Campos, T. Economon, T. Lukaczyk, S. Padron +# \version 7.0.4 "Blackbird" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md) +# +# SU2 is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# SU2 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with SU2. If not, see . + +# make print(*args) function available in PY2.6+, does'nt work on PY < 2.6 +from __future__ import print_function + +import sys +from TestCase import TestCase + +def main(): + '''This program runs SU2 and ensures that the output matches specified values. + This will be used to do checks when code is pushed to github + to make sure nothing is broken. ''' + + test_list = [] + + ########################## + ### Compressible Euler ### + ########################## + + # Channel + channel = TestCase('channel') + channel.cfg_dir = "euler/channel" + channel.cfg_file = "inv_channel_RK.cfg" + channel.test_iter = 20 + channel.test_vals = [-2.667326, 2.797439, 0.018717, 0.006906] + test_list.append(channel) + + # NACA0012 + naca0012 = TestCase('naca0012') + naca0012.cfg_dir = "euler/naca0012" + naca0012.cfg_file = "inv_NACA0012_Roe.cfg" + naca0012.test_iter = 20 + naca0012.test_vals = [-4.023999, -3.515034, 0.339426, 0.022217] + test_list.append(naca0012) + + # Supersonic wedge + wedge = TestCase('wedge') + wedge.cfg_dir = "euler/wedge" + wedge.cfg_file = "inv_wedge_HLLC.cfg" + wedge.test_iter = 20 + wedge.test_vals = [-0.942862, 4.784581, -0.208106, 0.036665] + test_list.append(wedge) + + # ONERA M6 Wing + oneram6 = TestCase('oneram6') + oneram6.cfg_dir = "euler/oneram6" + oneram6.cfg_file = "inv_ONERAM6.cfg" + oneram6.test_iter = 10 + oneram6.test_vals = [0.281704, 0.011821] + test_list.append(oneram6) + + # Fixed CL NACA0012 + fixedCL_naca0012 = TestCase('fixedcl_naca0012') + fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" + fixedCL_naca0012.cfg_file = "inv_NACA0012.cfg" + fixedCL_naca0012.test_iter = 10 + fixedCL_naca0012.test_vals = [-12.130189, -6.702728, 0.300000, 0.019470] + test_list.append(fixedCL_naca0012) + + # HYPERSONIC FLOW PAST BLUNT BODY + bluntbody = TestCase('bluntbody') + bluntbody.cfg_dir = "euler/bluntbody" + bluntbody.cfg_file = "blunt.cfg" + bluntbody.test_iter = 20 + bluntbody.test_vals = [0.540010, 6.916656, 0.000027, 1.869004] + test_list.append(bluntbody) + + ########################## + ### Compressible N-S ### + ########################## + + # Laminar flat plate + flatplate = TestCase('flatplate') + flatplate.cfg_dir = "navierstokes/flatplate" + flatplate.cfg_file = "lam_flatplate.cfg" + flatplate.test_iter = 20 + flatplate.test_vals = [-4.680775, 0.781235, -0.135957, 0.022978] + test_list.append(flatplate) + + # Laminar cylinder (steady) + cylinder = TestCase('cylinder') + cylinder.cfg_dir = "navierstokes/cylinder" + cylinder.cfg_file = "lam_cylinder.cfg" + cylinder.test_iter = 25 + cylinder.test_vals = [-6.765432, -1.297428, 0.019596, 0.310240] + test_list.append(cylinder) + + # Laminar cylinder (low Mach correction) + cylinder_lowmach = TestCase('cylinder_lowmach') + cylinder_lowmach.cfg_dir = "navierstokes/cylinder" + cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg" + cylinder_lowmach.test_iter = 25 + cylinder_lowmach.test_vals = [-6.850130, -1.388096, -0.056203, 108.140819] + test_list.append(cylinder_lowmach) + + # 2D Poiseuille flow (body force driven with periodic inlet / outlet) + poiseuille = TestCase('poiseuille') + poiseuille.cfg_dir = "navierstokes/poiseuille" + poiseuille.cfg_file = "lam_poiseuille.cfg" + poiseuille.test_iter = 10 + poiseuille.test_vals = [-5.048282, 0.650814, 0.008714, 13.677678] + test_list.append(poiseuille) + + # 2D Poiseuille flow (inlet profile file) + poiseuille_profile = TestCase('poiseuille_profile') + poiseuille_profile.cfg_dir = "navierstokes/poiseuille" + poiseuille_profile.cfg_file = "profile_poiseuille.cfg" + poiseuille_profile.test_iter = 10 + poiseuille_profile.test_vals = [-12.494721, -7.712408, -0.000000, 2.085796] + test_list.append(poiseuille_profile) + + ########################## + ### Compressible RANS ### + ########################## + + # RAE2822 SA + rae2822_sa = TestCase('rae2822_sa') + rae2822_sa.cfg_dir = "rans/rae2822" + rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg" + rae2822_sa.test_iter = 20 + rae2822_sa.test_vals = [-2.021218, -5.268447, 0.807465, 0.060897] + test_list.append(rae2822_sa) + + # RAE2822 SST + rae2822_sst = TestCase('rae2822_sst') + rae2822_sst.cfg_dir = "rans/rae2822" + rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" + rae2822_sst.test_iter = 20 + rae2822_sst.test_vals = [-0.510637, 4.876603, 0.812485, 0.061969] + test_list.append(rae2822_sst) + + # RAE2822 SST_SUST + rae2822_sst_sust = TestCase('rae2822_sst_sust') + rae2822_sst_sust.cfg_dir = "rans/rae2822" + rae2822_sst_sust.cfg_file = "turb_SST_SUST_RAE2822.cfg" + rae2822_sst_sust.test_iter = 20 + rae2822_sst_sust.test_vals = [-2.429813, 4.876602, 0.812485, 0.061969] + test_list.append(rae2822_sst_sust) + + # Flat plate + turb_flatplate = TestCase('turb_flatplate') + turb_flatplate.cfg_dir = "rans/flatplate" + turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" + turb_flatplate.test_iter = 20 + turb_flatplate.test_vals = [-4.157169, -6.737133, -0.176253, 0.057446] + test_list.append(turb_flatplate) + + # ONERA M6 Wing + turb_oneram6 = TestCase('turb_oneram6') + turb_oneram6.cfg_dir = "rans/oneram6" + turb_oneram6.cfg_file = "turb_ONERAM6.cfg" + turb_oneram6.test_iter = 10 + turb_oneram6.test_vals = [-2.372346, -6.579370, 0.229866, 0.147638] + test_list.append(turb_oneram6) + + # NACA0012 (SA, FUN3D finest grid results: CL=1.0983, CD=0.01242) + turb_naca0012_sa = TestCase('turb_naca0012_sa') + turb_naca0012_sa.cfg_dir = "rans/naca0012" + turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" + turb_naca0012_sa.test_iter = 10 + turb_naca0012_sa.test_vals = [-12.076819, -16.049252, 1.064326, 0.019770] + test_list.append(turb_naca0012_sa) + + # NACA0012 (SST, FUN3D finest grid results: CL=1.0840, CD=0.01253) + turb_naca0012_sst = TestCase('turb_naca0012_sst') + turb_naca0012_sst.cfg_dir = "rans/naca0012" + turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" + turb_naca0012_sst.test_iter = 10 + turb_naca0012_sst.test_vals = [-15.273728, -6.243783, 1.049988, 0.019165] + test_list.append(turb_naca0012_sst) + + # NACA0012 (SST_SUST, FUN3D finest grid results: CL=1.0840, CD=0.01253) + turb_naca0012_sst_sust = TestCase('turb_naca0012_sst_sust') + turb_naca0012_sst_sust.cfg_dir = "rans/naca0012" + turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg" + turb_naca0012_sst_sust.test_iter = 10 + turb_naca0012_sst_sust.test_vals = [-14.851214, -6.062566, 1.005233, 0.019014] + test_list.append(turb_naca0012_sst_sust) + + # PROPELLER + propeller = TestCase('propeller') + propeller.cfg_dir = "rans/propeller" + propeller.cfg_file = "propeller.cfg" + propeller.test_iter = 10 + propeller.test_vals = [-3.389576, -8.409529, 0.000048, 0.056329] + test_list.append(propeller) + + ################################# + ## Compressible RANS Restart ### + ################################# + + # NACA0012 SST Multigrid restart + turb_naca0012_sst_restart_mg = TestCase('turb_naca0012_sst_restart_mg') + turb_naca0012_sst_restart_mg.cfg_dir = "rans/naca0012" + turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" + turb_naca0012_sst_restart_mg.test_iter = 20 + turb_naca0012_sst_restart_mg.ntest_vals = 5 + turb_naca0012_sst_restart_mg.test_vals = [-7.653296, -7.729472, -1.981061, -0.000016, 0.079062] + test_list.append(turb_naca0012_sst_restart_mg) + + ############################# + ### Compressibele RANS UQ ### + ############################# + + # NACA0012 1c + turb_naca0012_1c = TestCase('turb_naca0012_1c') + turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" + turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" + turb_naca0012_1c.test_iter = 10 + turb_naca0012_1c.test_vals = [-4.906563, 1.338083, 6.086180, 2.412997] + test_list.append(turb_naca0012_1c) + + # NACA0012 2c + turb_naca0012_2c = TestCase('turb_naca0012_2c') + turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" + turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" + turb_naca0012_2c.test_iter = 10 + turb_naca0012_2c.test_vals = [-5.230218, 1.262228, 6.086157, 2.412759] + test_list.append(turb_naca0012_2c) + + # NACA0012 3c + turb_naca0012_3c = TestCase('turb_naca0012_3c') + turb_naca0012_3c.cfg_dir = "rans_uq/naca0012" + turb_naca0012_3c.cfg_file = "turb_NACA0012_uq_3c.cfg" + turb_naca0012_3c.test_iter = 10 + turb_naca0012_3c.test_vals = [-5.277130, 1.246265, 6.086053, 2.412464] + test_list.append(turb_naca0012_3c) + + # NACA0012 p1c1 + turb_naca0012_p1c1 = TestCase('turb_naca0012_p1c1') + turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" + turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" + turb_naca0012_p1c1.test_iter = 10 + turb_naca0012_p1c1.test_vals = [-5.012946, 1.309725, 6.086066, 2.412931] + test_list.append(turb_naca0012_p1c1) + + # NACA0012 p1c2 + turb_naca0012_p1c2 = TestCase('turb_naca0012_p1c2') + turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" + turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" + turb_naca0012_p1c2.test_iter = 10 + turb_naca0012_p1c2.test_vals = [-5.264019, 1.251314, 6.086451, 2.413070] + test_list.append(turb_naca0012_p1c2) + + ###################################### + ### Harmonic Balance ### + ###################################### + + # Description of the regression test + harmonic_balance = TestCase('harmonic_balance') + harmonic_balance.cfg_dir = "harmonic_balance" + harmonic_balance.cfg_file = "HB.cfg" + harmonic_balance.test_iter = 25 + harmonic_balance.test_vals = [-1.589863, 3.922098, -0.001443, 0.099457] + harmonic_balance.new_output = False + test_list.append(harmonic_balance) + + # Turbulent pitching NACA 64a010 airfoil + hb_rans_preconditioning = TestCase('hb_rans_preconditioning') + hb_rans_preconditioning.cfg_dir = "harmonic_balance/hb_rans_preconditioning" + hb_rans_preconditioning.cfg_file = "davis.cfg" + hb_rans_preconditioning.test_iter = 25 + hb_rans_preconditioning.test_vals = [-1.909633, -5.954752, 0.007773, 0.131217] + hb_rans_preconditioning.new_output = False + test_list.append(hb_rans_preconditioning) + + ###################################### + ### Moving Wall ### + ###################################### + + # Lid-driven cavity + cavity = TestCase('cavity') + cavity.cfg_dir = "moving_wall/cavity" + cavity.cfg_file = "lam_cavity.cfg" + cavity.test_iter = 25 + cavity.test_vals = [-5.627934, -0.164469, 0.051998, 2.547065] + test_list.append(cavity) + + # Spinning cylinder + spinning_cylinder = TestCase('spinning_cylinder') + spinning_cylinder.cfg_dir = "moving_wall/spinning_cylinder" + spinning_cylinder.cfg_file = "spinning_cylinder.cfg" + spinning_cylinder.test_iter = 25 + spinning_cylinder.test_vals = [-7.996313, -2.601764, 1.510692, 1.493876] + test_list.append(spinning_cylinder) + + ###################################### + ### Unsteady ### + ###################################### + + # Square cylinder + square_cylinder = TestCase('square_cylinder') + square_cylinder.cfg_dir = "unsteady/square_cylinder" + square_cylinder.cfg_file = "turb_square.cfg" + square_cylinder.test_iter = 3 + square_cylinder.test_vals = [-1.162572, 0.066371, 1.399790, 2.220393] + square_cylinder.unsteady = True + test_list.append(square_cylinder) + + # Gust + sine_gust = TestCase('sine_gust') + sine_gust.cfg_dir = "gust" + sine_gust.cfg_file = "inv_gust_NACA0012.cfg" + sine_gust.test_iter = 5 + sine_gust.test_vals = [-1.977545, 3.481778, -0.001667, -0.007429] + sine_gust.unsteady = True + test_list.append(sine_gust) + + # Aeroelastic + aeroelastic = TestCase('aeroelastic') + aeroelastic.cfg_dir = "aeroelastic" + aeroelastic.cfg_file = "aeroelastic_NACA64A010.cfg" + aeroelastic.test_iter = 2 + aeroelastic.test_vals = [0.079009, 0.033187, -0.001665, -0.000156] + aeroelastic.unsteady = True + test_list.append(aeroelastic) + + # Delayed Detached Eddy Simulation + ddes_flatplate = TestCase('ddes_flatplate') + ddes_flatplate.cfg_dir = "ddes/flatplate" + ddes_flatplate.cfg_file = "ddes_flatplate.cfg" + ddes_flatplate.test_iter = 10 + ddes_flatplate.test_vals = [-2.714758, -5.883004, -0.215005, 0.023783] + ddes_flatplate.unsteady = True + test_list.append(ddes_flatplate) + + ###################################### + ### NICFD ### + ###################################### + + # Rarefaction shock wave edge_VW + edge_VW = TestCase('edge_VW') + edge_VW.cfg_dir = "nicf/edge" + edge_VW.cfg_file = "edge_VW.cfg" + edge_VW.test_iter = 100 + edge_VW.test_vals = [-5.040283, 1.124491, -0.000009, 0.000000] + test_list.append(edge_VW) + + # Rarefaction shock wave edge_PPR + edge_PPR = TestCase('edge_PPR') + edge_PPR.cfg_dir = "nicf/edge" + edge_PPR.cfg_file = "edge_PPR.cfg" + edge_PPR.test_iter = 100 + edge_PPR.test_vals = [-5.401640, 0.738165, -0.000035, 0.000000] + test_list.append(edge_PPR) + + ###################################### + ### Turbomachinery ### + ###################################### + + # Jones APU Turbocharger + Jones_tc = TestCase('jones_turbocharger') + Jones_tc.cfg_dir = "turbomachinery/APU_turbocharger" + Jones_tc.cfg_file = "Jones.cfg" + Jones_tc.test_iter = 5 + Jones_tc.test_vals = [-5.280316, 0.379651, 44.725470, 2.271540] + Jones_tc.new_output = False + test_list.append(Jones_tc) + + # Jones APU Turbocharger restart + Jones_tc_rst = TestCase('jones_turbocharger_restart') + Jones_tc_rst.cfg_dir = "turbomachinery/APU_turbocharger" + Jones_tc_rst.cfg_file = "Jones_rst.cfg" + Jones_tc_rst.test_iter = 5 + Jones_tc_rst.test_vals = [-4.625319, -1.569634, 34.014100, 10.187660] + Jones_tc_rst.new_output = False + test_list.append(Jones_tc_rst) + + # 2D axial stage + axial_stage2D = TestCase('axial_stage2D') + axial_stage2D.cfg_dir = "turbomachinery/axial_stage_2D" + axial_stage2D.cfg_file = "Axial_stage2D.cfg" + axial_stage2D.test_iter = 20 + axial_stage2D.test_vals = [-1.933199, 5.381560, 73.357900, 1.780500] + axial_stage2D.new_output = False + test_list.append(axial_stage2D) + + # 2D transonic stator + transonic_stator = TestCase('transonic_stator') + transonic_stator.cfg_dir = "turbomachinery/transonic_stator_2D" + transonic_stator.cfg_file = "transonic_stator.cfg" + transonic_stator.test_iter = 20 + transonic_stator.test_vals = [-0.563540, 5.823232, 96.736080, 0.062426] + transonic_stator.new_output = False + test_list.append(transonic_stator) + + # 2D transonic stator restart + transonic_stator_rst = TestCase('transonic_stator_restart') + transonic_stator_rst.cfg_dir = "turbomachinery/transonic_stator_2D" + transonic_stator_rst.cfg_file = "transonic_stator_rst.cfg" + transonic_stator_rst.test_iter = 20 + transonic_stator_rst.test_vals = [-6.621626, -0.614379, 5.002986, 0.002951] + transonic_stator_rst.new_output = False + test_list.append(transonic_stator_rst) + + ###################################### + ### Sliding Mesh ### + ###################################### + + # Uniform flow + uniform_flow = TestCase('uniform_flow') + uniform_flow.cfg_dir = "sliding_interface/uniform_flow" + uniform_flow.cfg_file = "uniform_NN.cfg" + uniform_flow.test_iter = 5 + uniform_flow.test_vals = [5.000000, 0.000000, -0.188748, -10.631530] + uniform_flow.unsteady = True + uniform_flow.multizone = True + test_list.append(uniform_flow) + + # Channel_2D + channel_2D = TestCase('channel_2D') + channel_2D.cfg_dir = "sliding_interface/channel_2D" + channel_2D.cfg_file = "channel_2D_WA.cfg" + channel_2D.test_iter = 2 + channel_2D.test_vals = [2.000000, 0.000000, 0.398089, 0.352762, 0.405397] + channel_2D.unsteady = True + channel_2D.multizone = True + test_list.append(channel_2D) + + # Channel_3D + channel_3D = TestCase('channel_3D') + channel_3D.cfg_dir = "sliding_interface/channel_3D" + channel_3D.cfg_file = "channel_3D_WA.cfg" + channel_3D.test_iter = 2 + channel_3D.test_vals = [2.000000, 0.000000, 0.620151, 0.505157, 0.415249] + channel_3D.unsteady = True + channel_3D.multizone = True + test_list.append(channel_3D) + + # Pipe + pipe = TestCase('pipe') + pipe.cfg_dir = "sliding_interface/pipe" + pipe.cfg_file = "pipe_NN.cfg" + pipe.test_iter = 2 + pipe.test_vals = [0.150024, 0.491949, 0.677757, 0.963990, 1.006944] + pipe.unsteady = True + pipe.multizone = True + test_list.append(pipe) + + # Rotating cylinders + rotating_cylinders = TestCase('rotating_cylinders') + 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 = [3.000000, 0.000000, 0.777267, 1.134746, 1.224125] + rotating_cylinders.unsteady = True + rotating_cylinders.multizone = True + test_list.append(rotating_cylinders) + + # Supersonic vortex shedding + supersonic_vortex_shedding = TestCase('supersonic_vortex_shedding') + 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 = [5.000000, 0.000000, 1.216554, 1.639121] + supersonic_vortex_shedding.unsteady = True + supersonic_vortex_shedding.multizone = True + test_list.append(supersonic_vortex_shedding) + + # Bars_SST_2D + bars_SST_2D = TestCase('bars_SST_2D') + bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" + bars_SST_2D.cfg_file = "bars.cfg" + bars_SST_2D.test_iter = 13 + bars_SST_2D.test_vals = [13.000000, -0.619179, -1.564701] + bars_SST_2D.multizone = True + test_list.append(bars_SST_2D) + + ########################## + ### FEA - FSI ### + ########################## + + # Static beam, 3d + statbeam3d = TestCase('statbeam3d') + statbeam3d.cfg_dir = "fea_fsi/StatBeam_3d" + statbeam3d.cfg_file = "configBeam_3d.cfg" + statbeam3d.test_iter = 0 + statbeam3d.test_vals = [-2.378370, -1.585252, -2.028505, 6.4359e+04] + test_list.append(statbeam3d) + + # Dynamic beam, 2d + dynbeam2d = TestCase('dynbeam2d') + dynbeam2d.cfg_dir = "fea_fsi/DynBeam_2d" + dynbeam2d.cfg_file = "configBeam_2d.cfg" + dynbeam2d.test_iter = 6 + dynbeam2d.test_vals = [-3.240015, 2.895057, -0.353146, 66127.000000] + dynbeam2d.unsteady = True + test_list.append(dynbeam2d) + + # FSI, 2d + fsi2d = TestCase('fsi2d') + fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" + fsi2d.cfg_file = "configFSI.cfg" + fsi2d.test_iter = 4 + fsi2d.test_vals = [4.000000, 0.000000, -3.764089, -4.081119] + fsi2d.multizone= True + fsi2d.unsteady = True + test_list.append(fsi2d) + + # FSI, Static, 2D, new mesh solver + stat_fsi = TestCase('stat_fsi') + stat_fsi.cfg_dir = "fea_fsi/stat_fsi" + stat_fsi.cfg_file = "config.cfg" + stat_fsi.test_iter = 7 + stat_fsi.test_vals = [-3.242834, -4.866608, 0.000000, 11.000000] + stat_fsi.multizone = True + test_list.append(stat_fsi) + + # FSI, Dynamic, 2D, new mesh solver + dyn_fsi = TestCase('dyn_fsi') + dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" + dyn_fsi.cfg_file = "config.cfg" + dyn_fsi.test_iter = 4 + dyn_fsi.test_vals = [-4.379823, -4.005990, 0.000000, 74.000000] + dyn_fsi.multizone = True + dyn_fsi.unsteady = True + test_list.append(dyn_fsi) + + # FSI, Static, 2D, new mesh solver, restart + stat_fsi_restart = TestCase('stat_fsi_restart') + stat_fsi_restart.cfg_dir = "fea_fsi/stat_fsi" + stat_fsi_restart.cfg_file = "config_restart.cfg" + stat_fsi_restart.test_iter = 1 + stat_fsi_restart.test_vals = [-3.474239, -4.250710, 0.000000, 36.000000] + stat_fsi_restart.multizone = True + test_list.append(stat_fsi_restart) + + ############################################## + ### Method of Manufactured Solutions (MMS) ### + ############################################## + + # FVM, compressible, laminar N-S + mms_fvm_ns = TestCase('mms_fvm_ns') + mms_fvm_ns.cfg_dir = "mms/fvm_navierstokes" + mms_fvm_ns.cfg_file = "lam_mms_roe.cfg" + mms_fvm_ns.test_iter = 20 + mms_fvm_ns.test_vals = [-2.851428, 2.192348, 0.000000, 0.000000] + test_list.append(mms_fvm_ns) + + ###################################### + ### RUN TESTS ### + ###################################### + + for test in test_list: + test.su2_exec = "SU2_CFD -t 2" + test.timeout = 600 + test.tol = 1e-4 + #end + + pass_list = [ test.run_test() for test in test_list ] + + # Tests summary + print('==================================================================') + print('Summary of the parallel tests') + print('python version:', sys.version) + for i, test in enumerate(test_list): + if (pass_list[i]): + print(' passed - %s'%test.tag) + else: + print('* FAILED - %s'%test.tag) + + if all(pass_list): + sys.exit(0) + else: + sys.exit(1) + # done + +if __name__ == '__main__': + main() diff --git a/UnitTests/Common/simple_ad_test.cpp b/UnitTests/Common/simple_ad_test.cpp index 96777f1913ce..dc6075b3028a 100644 --- a/UnitTests/Common/simple_ad_test.cpp +++ b/UnitTests/Common/simple_ad_test.cpp @@ -29,7 +29,7 @@ #include "catch.hpp" -#include "../../Common/include/datatype_structure.hpp" +#include "../../Common/include/basic_types/datatype_structure.hpp" su2double func(const su2double& x) { return x * x * x; diff --git a/UnitTests/Common/simple_directdiff_test.cpp b/UnitTests/Common/simple_directdiff_test.cpp index 04bfc18b0d66..48007f7e1e27 100644 --- a/UnitTests/Common/simple_directdiff_test.cpp +++ b/UnitTests/Common/simple_directdiff_test.cpp @@ -29,7 +29,7 @@ #include "catch.hpp" -#include "../../Common/include/datatype_structure.hpp" +#include "../../Common/include/basic_types/datatype_structure.hpp" su2double func(const su2double& x) { return x * x * x;