diff --git a/Common/include/omp_structure.hpp b/Common/include/omp_structure.hpp index 25b8e49d2a02..9a610d49f26e 100644 --- a/Common/include/omp_structure.hpp +++ b/Common/include/omp_structure.hpp @@ -63,6 +63,11 @@ */ inline constexpr int omp_get_max_threads(void) {return 1;} +/*! + * \brief Number of threads in current team. + */ +inline constexpr int omp_get_num_threads(void) {return 1;} + /*! * \brief Index of current thread, akin to MPI rank. */ diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 759b64129317..6528ed174121 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1933,6 +1933,7 @@ static const map Input_Ref_Map = CCreateMap +struct C2DDummyLastView +{ + using Index = typename T::Index; + using Scalar = typename T::Scalar; + + T& data; + + C2DDummyLastView() = delete; + + C2DDummyLastView(T& ref) : data(ref) {} + + template::value, bool>::type = 0> + Scalar& operator() (Index i, Index) noexcept + { + return data(i); + } + + const Scalar& operator() (Index i, Index) const noexcept + { + return data(i); + } +}; + +/*! + * \class C3DDummyMiddleView + * \brief Helper class, adds dummy middle dimension to a reference of a + * matrix object making it a dummy 3D array. + * \note The constness of the object is derived from the template type, but + * we allways keep a reference, never a copy of the associated matrix. + */ +template +struct C3DDummyMiddleView +{ + using Index = typename T::Index; + using Scalar = typename T::Scalar; + + T& data; + + C3DDummyMiddleView() = delete; + + C3DDummyMiddleView(T& ref) : data(ref) {} + + template::value, bool>::type = 0> + Scalar& operator() (Index i, Index, Index k) noexcept + { + return data(i,k); + } + + const Scalar& operator() (Index i, Index, Index k) const noexcept + { + return data(i,k); + } +}; + /*! * \brief Useful typedefs with default template parameters */ diff --git a/Common/src/meson.build b/Common/src/meson.build index ea7098f10a8d..9aabdeb41a80 100644 --- a/Common/src/meson.build +++ b/Common/src/meson.build @@ -35,8 +35,9 @@ if get_option('enable-normal') common = static_library('SU2Common', common_src, install : false, - dependencies : su2_deps, - cpp_args: [default_warning_flags, su2_cpp_args]) + dependencies : su2_deps, + cpp_args: [default_warning_flags, su2_cpp_args]) + common_dep = declare_dependency(link_with: common, include_directories : common_include) endif @@ -46,8 +47,8 @@ if get_option('enable-autodiff') commonAD = static_library('SU2CommonAD', common_src, install : false, - dependencies : [su2_deps, codi_dep], - cpp_args: [default_warning_flags, su2_cpp_args, codi_rev_args]) + dependencies : [su2_deps, codi_dep], + cpp_args: [default_warning_flags, su2_cpp_args, codi_rev_args]) commonAD_dep = declare_dependency(link_with: commonAD, include_directories : common_include) @@ -59,8 +60,8 @@ if get_option('enable-directdiff') commonDD = static_library('SU2CommonDD', common_src, install : false, - dependencies : [su2_deps, codi_dep], - cpp_args: [default_warning_flags, su2_cpp_args, codi_for_args]) + dependencies : [su2_deps, codi_dep], + cpp_args: [default_warning_flags, su2_cpp_args, codi_for_args]) commonDD_dep = declare_dependency(link_with: commonDD, include_directories : common_include) diff --git a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp new file mode 100644 index 000000000000..98d152ef8dd5 --- /dev/null +++ b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp @@ -0,0 +1,189 @@ +/*! + * \file computeGradientsGreenGauss.hpp + * \brief Generic implementation of Green-Gauss gradient computation. + * \note This allows the same implementation to be used for conservative + * and primitive variables of any solver. + * \author P. Gomes + * \version 7.0.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2019, 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 . + */ + +#include "../../../Common/include/omp_structure.hpp" + + +/*! + * \brief Compute the gradient of a field using the Green-Gauss theorem. + * \note Gradients can be computed only for a contiguous range of variables, defined + * by [varBegin, varEnd[ (e.g. 0,1 computes the gradient of the 1st variable). + * This can be used, for example, to compute only velocity gradients. + * \note The function uses an optional solver object to perform communications, if + * none (nullptr) is provided the function does not fail (the objective of + * this is to improve test-ability). + * \param[in] solver - Optional, solver associated with the field (used only for MPI). + * \param[in] kindMpiComm - Type of MPI communication required. + * \param[in] kindPeriodicComm - Type of periodic communication required. + * \param[in] geometry - Geometric grid properties. + * \param[in] config - Configuration of the problem, used to identify types of boundaries. + * \param[in] field - Generic object implementing operator (iPoint, iVar). + * \param[in] varBegin - Index of first variable for which to compute the gradient. + * \param[in] varEnd - Index of last variable for which to compute the gradient. + * \param[out] gradient - Generic object implementing operator (iPoint, iVar, iDim). + */ +template +void computeGradientsGreenGauss(CSolver* solver, + MPI_QUANTITIES kindMpiComm, + PERIODIC_QUANTITIES kindPeriodicComm, + CGeometry& geometry, + CConfig& config, + const FieldType& field, + size_t varBegin, + size_t varEnd, + GradientType& gradient) +{ + size_t nPointDomain = geometry.GetnPointDomain(); + size_t nDim = geometry.GetnDim(); + +#ifdef HAVE_OMP + constexpr size_t OMP_MAX_CHUNK = 512; + + size_t chunkSize = computeStaticChunkSize(nPointDomain, + omp_get_max_threads(), OMP_MAX_CHUNK); +#endif + + /*--- Start OpenMP parallel section. ---*/ + + SU2_OMP_PARALLEL + { + /*--- For each (non-halo) volume integrate over its faces (edges). ---*/ + + SU2_OMP_FOR_DYN(chunkSize) + for (size_t iPoint = 0; iPoint < nPointDomain; ++iPoint) + { + auto node = geometry.node[iPoint]; + + AD::StartPreacc(); + AD::SetPreaccIn(node->GetVolume()); + AD::SetPreaccIn(node->GetPeriodicVolume()); + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + AD::SetPreaccIn(field(iPoint,iVar)); + + /*--- Clear the gradient. --*/ + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + for (size_t iDim = 0; iDim < nDim; ++iDim) + gradient(iPoint, iVar, iDim) = 0.0; + + /*--- Handle averaging and division by volume in one constant. ---*/ + + su2double halfOnVol = 0.5 / (node->GetVolume()+node->GetPeriodicVolume()); + + /*--- Add a contribution due to each neighbor. ---*/ + + for (size_t iNeigh = 0; iNeigh < node->GetnPoint(); ++iNeigh) + { + size_t iEdge = node->GetEdge(iNeigh); + size_t jPoint = node->GetPoint(iNeigh); + + /*--- Determine if edge points inwards or outwards of iPoint. + * If inwards we need to flip the area vector. ---*/ + + su2double dir = (iPoint == geometry.edge[iEdge]->GetNode(0))? 1.0 : -1.0; + su2double weight = dir * halfOnVol; + + const su2double* area = geometry.edge[iEdge]->GetNormal(); + AD::SetPreaccIn(area, nDim); + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + { + AD::SetPreaccIn(field(jPoint,iVar)); + + su2double flux = weight * (field(iPoint,iVar) + field(jPoint,iVar)); + + for (size_t iDim = 0; iDim < nDim; ++iDim) + gradient(iPoint, iVar, iDim) += flux * area[iDim]; + } + + } + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + for (size_t iDim = 0; iDim < nDim; ++iDim) + AD::SetPreaccOut(gradient(iPoint,iVar,iDim)); + + AD::EndPreacc(); + } + + /*--- Add boundary fluxes. ---*/ + + for (size_t iMarker = 0; iMarker < geometry.GetnMarker(); ++iMarker) + { + if ((config.GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && + (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) + { + /*--- Work is shared in inner loop as two markers + * may try to update the same point. ---*/ + + SU2_OMP_FOR_STAT(OMP_MAX_CHUNK) + for (size_t iVertex = 0; iVertex < geometry.GetnVertex(iMarker); ++iVertex) + { + size_t iPoint = geometry.vertex[iMarker][iVertex]->GetNode(); + auto node = geometry.node[iPoint]; + + /*--- Halo points do not need to be considered. ---*/ + + if (!node->GetDomain()) continue; + + su2double volume = node->GetVolume() + node->GetPeriodicVolume(); + + const su2double* area = geometry.vertex[iMarker][iVertex]->GetNormal(); + + for (size_t iVar = varBegin; iVar < varEnd; iVar++) + { + su2double flux = field(iPoint,iVar) / volume; + + for (size_t iDim = 0; iDim < nDim; iDim++) + gradient(iPoint, iVar, iDim) -= flux * area[iDim]; + } + } + } + } + + } // end SU2_OMP_PARALLEL + + /*--- If no solver was provided we do not communicate ---*/ + + if (solver == nullptr) return; + + /*--- Account for periodic contributions. ---*/ + + for (size_t iPeriodic = 1; iPeriodic <= config.GetnMarker_Periodic()/2; ++iPeriodic) + { + solver->InitiatePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm); + solver->CompletePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm); + } + + /*--- Obtain the gradients at halo points from the MPI ranks that own them. ---*/ + + solver->InitiateComms(&geometry, &config, kindMpiComm); + solver->CompleteComms(&geometry, &config, kindMpiComm); + +} diff --git a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp new file mode 100644 index 000000000000..1f1766813662 --- /dev/null +++ b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp @@ -0,0 +1,315 @@ +/*! + * \file computeGradientsLeastSquares.hpp + * \brief Generic implementation of Least-Squares gradient computation. + * \note This allows the same implementation to be used for conservative + * and primitive variables of any solver. + * \version 7.0.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2019, 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 . + */ + +#include "../../../Common/include/omp_structure.hpp" + + +/*! + * \brief Compute the gradient of a field using inverse-distance-weighted or + * unweighted Least-Squares approximation. + * \note See notes from computeGradientsGreenGauss.hpp. + * \param[in] solver - Optional, solver associated with the field (used only for MPI). + * \param[in] kindMpiComm - Type of MPI communication required. + * \param[in] kindPeriodicComm - Type of periodic communication required. + * \param[in] geometry - Geometric grid properties. + * \param[in] weighted - Use inverse-distance weights. + * \param[in] config - Configuration of the problem, used to identify types of boundaries. + * \param[in] field - Generic object implementing operator (iPoint, iVar). + * \param[in] varBegin - Index of first variable for which to compute the gradient. + * \param[in] varEnd - Index of last variable for which to compute the gradient. + * \param[out] gradient - Generic object implementing operator (iPoint, iVar, iDim). + * \param[out] Rmatrix - Generic object implementing operator (iPoint, iDim, iDim). + */ +template +void computeGradientsLeastSquares(CSolver* solver, + MPI_QUANTITIES kindMpiComm, + PERIODIC_QUANTITIES kindPeriodicComm, + CGeometry& geometry, + CConfig& config, + bool weighted, + const FieldType& field, + size_t varBegin, + size_t varEnd, + GradientType& gradient, + RMatrixType& Rmatrix) +{ + constexpr size_t MAXNDIM = 3; + + size_t nPointDomain = geometry.GetnPointDomain(); + size_t nDim = geometry.GetnDim(); + +#ifdef HAVE_OMP + constexpr size_t OMP_MAX_CHUNK = 512; + + size_t chunkSize = computeStaticChunkSize(nPointDomain, + omp_get_max_threads(), OMP_MAX_CHUNK); +#endif + + /*--- Start OpenMP parallel section. ---*/ + + SU2_OMP_PARALLEL + { + /*--- First loop over non-halo points of the grid. ---*/ + + SU2_OMP_FOR_DYN(chunkSize) + for (size_t iPoint = 0; iPoint < nPointDomain; ++iPoint) + { + auto node = geometry.node[iPoint]; + const su2double* coord_i = node->GetCoord(); + + AD::StartPreacc(); + AD::SetPreaccIn(coord_i, nDim); + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + AD::SetPreaccIn(field(iPoint,iVar)); + + /*--- Clear gradient and Rmatrix. ---*/ + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + for (size_t iDim = 0; iDim < nDim; ++iDim) + gradient(iPoint, iVar, iDim) = 0.0; + + for (size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t jDim = 0; jDim < nDim; ++jDim) + Rmatrix(iPoint, iDim, jDim) = 0.0; + + + for (size_t iNeigh = 0; iNeigh < node->GetnPoint(); ++iNeigh) + { + size_t jPoint = node->GetPoint(iNeigh); + + const su2double* coord_j = geometry.node[jPoint]->GetCoord(); + AD::SetPreaccIn(coord_j, nDim); + + /*--- Distance vector from iPoint to jPoint ---*/ + + su2double dist_ij[MAXNDIM]; + + for (size_t iDim = 0; iDim < nDim; ++iDim) + dist_ij[iDim] = coord_j[iDim] - coord_i[iDim]; + + /*--- Compute inverse weight, default 1 (unweighted). ---*/ + + su2double weight = 1.0; + + if (weighted) + { + weight = 0.0; + for (size_t iDim = 0; iDim < nDim; ++iDim) + weight += dist_ij[iDim] * dist_ij[iDim]; + } + + /*--- Sumations for entries of upper triangular matrix R. ---*/ + + if (weight > 0.0) + { + weight = 1.0 / weight; + + Rmatrix(iPoint,0,0) += dist_ij[0]*dist_ij[0]*weight; + Rmatrix(iPoint,0,1) += dist_ij[0]*dist_ij[1]*weight; + Rmatrix(iPoint,1,1) += dist_ij[1]*dist_ij[1]*weight; + + if (nDim == 3) + { + Rmatrix(iPoint,0,2) += dist_ij[0]*dist_ij[2]*weight; + Rmatrix(iPoint,1,2) += dist_ij[1]*dist_ij[2]*weight; + Rmatrix(iPoint,2,1) += dist_ij[0]*dist_ij[2]*weight; + Rmatrix(iPoint,2,2) += dist_ij[2]*dist_ij[2]*weight; + } + + /*--- Entries of c:= transpose(A)*b ---*/ + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + { + AD::SetPreaccIn(field(jPoint,iVar)); + + su2double delta_ij = weight * (field(jPoint,iVar) - field(iPoint,iVar)); + + for (size_t iDim = 0; iDim < nDim; ++iDim) + gradient(iPoint, iVar, iDim) += dist_ij[iDim] * delta_ij; + } + } + } + + for (size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t jDim = 0; jDim < nDim; ++jDim) + AD::SetPreaccOut(Rmatrix(iPoint, iDim, jDim)); + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + for (size_t iDim = 0; iDim < nDim; ++iDim) + AD::SetPreaccOut(gradient(iPoint, iVar, iDim)); + + AD::EndPreacc(); + } + + /*--- Correct the gradient values across any periodic boundaries. ---*/ + + if (solver != nullptr) + { + SU2_OMP_MASTER + { + for (size_t iPeriodic = 1; iPeriodic <= config.GetnMarker_Periodic()/2; ++iPeriodic) + { + solver->InitiatePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm); + solver->CompletePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm); + } + } + SU2_OMP_BARRIER + } + + /*--- Second loop over points of the grid to compute final gradient. ---*/ + + SU2_OMP_FOR_DYN(chunkSize) + for (size_t iPoint = 0; iPoint < nPointDomain; ++iPoint) + { + /*--- Entries of upper triangular matrix R. ---*/ + + su2double r11 = Rmatrix(iPoint,0,0); + su2double r12 = Rmatrix(iPoint,0,1); + su2double r22 = Rmatrix(iPoint,1,1); + su2double r13 = 0.0, r23 = 0.0, r23_a = 0.0, r23_b = 0.0, r33 = 0.0; + + AD::StartPreacc(); + AD::SetPreaccIn(r11); + AD::SetPreaccIn(r12); + AD::SetPreaccIn(r22); + + if (r11 >= 0.0) r11 = sqrt(r11); + if (r11 >= 0.0) r12 /= r11; else r12 = 0.0; + su2double tmp = r22-r12*r12; + if (tmp >= 0.0) r22 = sqrt(tmp); else r22 = 0.0; + + if (nDim == 3) { + r13 = Rmatrix(iPoint,0,2); + r23_a = Rmatrix(iPoint,1,2); + r23_b = Rmatrix(iPoint,2,1); + r33 = Rmatrix(iPoint,2,2); + + AD::SetPreaccIn(r13); + AD::SetPreaccIn(r23_a); + AD::SetPreaccIn(r23_b); + AD::SetPreaccIn(r33); + + if (r11 >= 0.0) r13 /= r11; else r13 = 0.0; + + if ((r22 >= 0.0) && (r11*r22 >= 0.0)) { + r23 = r23_a/r22 - r23_b*r12/(r11*r22); + } else { + r23 = 0.0; + } + + tmp = r33 - r23*r23 - r13*r13; + if (tmp >= 0.0) r33 = sqrt(tmp); else r33 = 0.0; + } + + /*--- Compute determinant ---*/ + + su2double detR2 = (r11*r22)*(r11*r22); + if (nDim == 3) detR2 *= r33*r33; + + /*--- Detect singular matrices ---*/ + + bool singular = false; + + if (detR2 <= EPS) { + detR2 = 1.0; + singular = true; + } + + /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ + + su2double Smatrix[MAXNDIM][MAXNDIM]; + + if (singular) { + for (size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t jDim = 0; jDim < nDim; ++jDim) + Smatrix[iDim][jDim] = 0.0; + } + else { + if (nDim == 2) { + Smatrix[0][0] = (r12*r12+r22*r22)/detR2; + Smatrix[0][1] = -r11*r12/detR2; + Smatrix[1][0] = Smatrix[0][1]; + Smatrix[1][1] = r11*r11/detR2; + } + else { + su2double z11 = r22*r33; + su2double z12 =-r12*r33; + su2double z13 = r12*r23-r13*r22; + su2double z22 = r11*r33; + su2double z23 =-r11*r23; + su2double z33 = r11*r22; + + Smatrix[0][0] = (z11*z11+z12*z12+z13*z13)/detR2; + Smatrix[0][1] = (z12*z22+z13*z23)/detR2; + Smatrix[0][2] = (z13*z33)/detR2; + Smatrix[1][0] = Smatrix[0][1]; + Smatrix[1][1] = (z22*z22+z23*z23)/detR2; + Smatrix[1][2] = (z23*z33)/detR2; + Smatrix[2][0] = Smatrix[0][2]; + Smatrix[2][1] = Smatrix[1][2]; + Smatrix[2][2] = (z33*z33)/detR2; + } + } + + for (size_t iDim = 0; iDim < nDim; ++iDim) + for (size_t jDim = 0; jDim < nDim; ++jDim) + AD::SetPreaccOut(Smatrix[iDim][jDim]); + + AD::EndPreacc(); + + /*--- Computation of the gradient: S*c ---*/ + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + { + su2double Cvector[MAXNDIM]; + + for (size_t iDim = 0; iDim < nDim; ++iDim) + { + Cvector[iDim] = 0.0; + for (size_t jDim = 0; jDim < nDim; ++jDim) + Cvector[iDim] += Smatrix[iDim][jDim] * gradient(iPoint, iVar, jDim); + } + + for (size_t iDim = 0; iDim < nDim; ++iDim) + gradient(iPoint, iVar, iDim) = Cvector[iDim]; + } + } + + } // end SU2_OMP_PARALLEL + + /*--- If no solver was provided we do not communicate ---*/ + + if (solver == nullptr) return; + + /*--- Obtain the gradients at halo points from the MPI ranks that own them. ---*/ + + solver->InitiateComms(&geometry, &config, kindMpiComm); + solver->CompleteComms(&geometry, &config, kindMpiComm); + +} diff --git a/SU2_CFD/include/limiters/CLimiterDetails.hpp b/SU2_CFD/include/limiters/CLimiterDetails.hpp new file mode 100644 index 000000000000..bd35109f1dd2 --- /dev/null +++ b/SU2_CFD/include/limiters/CLimiterDetails.hpp @@ -0,0 +1,335 @@ +/*! + * \file CLimiterDetails.hpp + * \brief A template class that allows defining limiters via + * specialization of particular details. + * \author P. Gomes + * \version 7.0.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2019, 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 . + */ + + +/*! + * \brief A traits class for limiters, see notes for "computeLimiters_impl()". + * \note There is no default implementation (the code will compile but not + * link) specialization is mandatory. + */ +template +struct CLimiterDetails +{ + /*! + * \brief Compute any global value that may be needed by the other functions. + * \note This function is called once by multiple threads. + */ + template + inline void preprocess(CGeometry&, CConfig&, size_t varBegin, + size_t varEnd, const FieldType&); + + /*! + * \brief Geometric modifier (e.g. increase limiting near sharp edges). + * \note This function is called once per point inside an AD pre- + * -accumulation region, newly used variables should be registered. + */ + inline su2double geometricFactor(size_t iPoint, CGeometry&) const; + + /*! + * \brief Smooth (usually) function of the maximum/minimum (positive/negative) + * gradient projections onto the edges, and the deltas over direct neighbors. + * Both proj and delta may be 0.0, beware of divisions. + * \note This function is called twice (min/max) per point per variable + * (also inside an AD pre-accumulation region). + */ + inline su2double limiterFunction(size_t iVar, su2double proj, su2double delta) const; +}; + + +/*! + * \brief Common small functions used by limiters. + */ +namespace LimiterHelpers +{ + inline passivedouble epsilon() {return std::numeric_limits::epsilon();} + + inline su2double venkatFunction(su2double proj, su2double delta, su2double eps2) + { + su2double y = delta*(delta+proj) + eps2; + return (y + delta*proj) / (y + 2*proj*proj); + } + + inline su2double raisedSine(su2double dist) + { + su2double factor = 0.5*(1.0+dist+sin(PI_NUMBER*dist)/PI_NUMBER); + return max(0.0, min(factor, 1.0)); + } +} + + +/*! + * \brief Barth-Jespersen specialization. + */ +template<> +struct CLimiterDetails +{ + su2double eps2; + + /*! + * \brief Set a small epsilon to avoid divisions by 0. + */ + template + inline void preprocess(Ts&...) {eps2 = LimiterHelpers::epsilon();} + + /*! + * \brief No geometric modification for this kind of limiter. + */ + template + inline su2double geometricFactor(Ts&...) const {return 1.0;} + + /*! + * \brief Venkatakrishnan function with a numerical epsilon. + */ + inline su2double limiterFunction(size_t, su2double proj, su2double delta) const + { + return LimiterHelpers::venkatFunction(proj, delta, eps2); + } +}; + + +/*! + * \brief Venkatakrishnan specialization. + */ +template<> +struct CLimiterDetails +{ + su2double eps2; + + /*! + * \brief Store the reference lenght based eps^2 parameter, + * limited to a small number to avoid divisions by 0. + */ + template + inline void preprocess(CGeometry&, CConfig& config, Ts&...) + { + su2double L = config.GetRefElemLength(); + su2double K = config.GetVenkat_LimiterCoeff(); + su2double eps1 = fabs(L*K); + eps2 = max(eps1*eps1*eps1, LimiterHelpers::epsilon()); + } + + /*! + * \brief No geometric modification for this kind of limiter. + */ + template + inline su2double geometricFactor(Ts&...) const {return 1.0;} + + /*! + * \brief Smooth function that disables limiting in smooth regions. + */ + inline su2double limiterFunction(size_t, su2double proj, su2double delta) const + { + return LimiterHelpers::venkatFunction(proj, delta, eps2); + } +}; + + +/*! + * \brief Venkatakrishnan-Wang specialization. + */ +template<> +struct CLimiterDetails +{ + su2activematrix fieldMin, fieldMax; + su2activevector eps2; + + /*! + * \brief Store the solution range based eps^2 parameter. + */ + template + inline void preprocess(CGeometry& geometry, CConfig& config, size_t varBegin, + size_t varEnd, const FieldType& field) + { + /*--- Determine the max and min global value for each variable. + * Each thread initially works on one row of fieldMin/Max (the + * rows are padded to a multiple of 8 to avoid false sharing), + * then the master thread performs a reduction over all threads + * and mpi ranks onto row 0, the final result. ---*/ + + size_t nThread = omp_get_num_threads(); + size_t nCols = roundUpDiv(varEnd,8)*8; + + SU2_OMP_MASTER + { + su2double largeNum = 0.1*std::numeric_limits::max(); + fieldMin.resize(nThread, nCols) = largeNum; + fieldMax.resize(nThread, nCols) =-largeNum; + eps2.resize(nCols); + } + SU2_OMP_BARRIER + + /*--- Per thread reduction. ---*/ + + SU2_OMP_FOR_STAT(512) + for(size_t iPoint = 0; iPoint < geometry.GetnPointDomain(); ++iPoint) + { + size_t iThread = omp_get_thread_num(); + + for(size_t iVar = varBegin; iVar < varEnd; ++iVar) + { + fieldMin(iThread, iVar) = min(fieldMin(iThread, iVar), field(iPoint, iVar)); + fieldMax(iThread, iVar) = max(fieldMax(iThread, iVar), field(iPoint, iVar)); + } + } + + SU2_OMP_MASTER + { + /*--- Per rank reduction. ---*/ + + for(size_t iThread = 1; iThread < nThread; ++iThread) + { + for(size_t iVar = varBegin; iVar < varEnd; ++iVar) + { + fieldMin(0,iVar) = min(fieldMin(0,iVar), fieldMin(iThread, iVar)); + fieldMax(0,iVar) = max(fieldMax(0,iVar), fieldMax(iThread, iVar)); + } + } + + /*--- Global reduction, (re)using eps2 as the recv buffer. ---*/ + + SU2_MPI::Allreduce(fieldMin[0], eps2.data(), nCols, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + + for(size_t iVar = varBegin; iVar < varEnd; ++iVar) + fieldMin(0,iVar) = eps2(iVar); + + SU2_MPI::Allreduce(fieldMax[0], eps2.data(), nCols, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + + for(size_t iVar = varBegin; iVar < varEnd; ++iVar) + fieldMax(0,iVar) = eps2(iVar); + + /*--- Compute eps^2 ---*/ + + su2double K = config.GetVenkat_LimiterCoeff(); + + for(size_t iVar = varBegin; iVar < varEnd; ++iVar) + { + su2double range = fieldMax(0,iVar) - fieldMin(0,iVar); + eps2(iVar) = max(pow(K*range, 2), LimiterHelpers::epsilon()); + } + } + SU2_OMP_BARRIER + + } + + /*! + * \brief No geometric modification for this kind of limiter. + */ + template + inline su2double geometricFactor(Ts&...) const {return 1.0;} + + /*! + * \brief Smooth function that disables limiting in smooth regions. + */ + inline su2double limiterFunction(size_t iVar, su2double proj, su2double delta) const + { + AD::SetPreaccIn(eps2(iVar)); + return LimiterHelpers::venkatFunction(proj, delta, eps2(iVar)); + } +}; + + +/*! + * \brief Venkatakrishnan with sharp edge modification. + */ +template<> +struct CLimiterDetails +{ + su2double eps1, eps2, sharpCoeff; + + /*! + * \brief Store the reference lenght based eps^2 parameter. + */ + template + inline void preprocess(CGeometry&, CConfig& config, Ts&...) + { + sharpCoeff = config.GetAdjSharp_LimiterCoeff(); + su2double L = config.GetRefElemLength(); + su2double K = config.GetVenkat_LimiterCoeff(); + eps1 = fabs(L*K); + eps2 = max(eps1*eps1*eps1, LimiterHelpers::epsilon()); + } + + /*! + * \brief Full limiting (1st order) near sharp edges. + */ + inline su2double geometricFactor(size_t iPoint, CGeometry& geometry) const + { + AD::SetPreaccIn(geometry.node[iPoint]->GetSharpEdge_Distance()); + su2double dist = geometry.node[iPoint]->GetSharpEdge_Distance()/(sharpCoeff*eps1)-1.0; + return LimiterHelpers::raisedSine(dist); + } + + /*! + * \brief Smooth function that disables limiting in smooth regions. + */ + inline su2double limiterFunction(size_t, su2double proj, su2double delta) const + { + return LimiterHelpers::venkatFunction(proj, delta, eps2); + } +}; + + +/*! + * \brief Venkatakrishnan with wall distance modification. + */ +template<> +struct CLimiterDetails +{ + su2double eps1, eps2, sharpCoeff; + + /*! + * \brief Store the reference lenght based eps^2 parameter. + */ + template + inline void preprocess(CGeometry&, CConfig& config, Ts&...) + { + sharpCoeff = config.GetAdjSharp_LimiterCoeff(); + su2double L = config.GetRefElemLength(); + su2double K = config.GetVenkat_LimiterCoeff(); + eps1 = fabs(L*K); + eps2 = max(eps1*eps1*eps1, LimiterHelpers::epsilon()); + } + + /*! + * \brief Full limiting (1st order) near walls. + */ + inline su2double geometricFactor(size_t iPoint, CGeometry& geometry) const + { + AD::SetPreaccIn(geometry.node[iPoint]->GetWall_Distance()); + su2double dist = geometry.node[iPoint]->GetWall_Distance()/(sharpCoeff*eps1)-1.0; + return LimiterHelpers::raisedSine(dist); + } + + /*! + * \brief Smooth function that disables limiting in smooth regions. + */ + inline su2double limiterFunction(size_t, su2double proj, su2double delta) const + { + return LimiterHelpers::venkatFunction(proj, delta, eps2); + } +}; diff --git a/SU2_CFD/include/limiters/computeLimiters.hpp b/SU2_CFD/include/limiters/computeLimiters.hpp new file mode 100644 index 000000000000..62b731467c77 --- /dev/null +++ b/SU2_CFD/include/limiters/computeLimiters.hpp @@ -0,0 +1,102 @@ +/*! + * \file computeLimiters.hpp + * \brief Compute limiters wrapper function. + * \author P. Gomes + * \version 7.0.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2019, 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 . + */ + +#include "CLimiterDetails.hpp" +#include "computeLimiters_impl.hpp" + +/*! + * \brief A wrapper funtion that calls specialized implementations depending + * on "LimiterKind". Those implementations are generated by instantiating + * versions of "computeLimiters_impl" with appropriate specializations + * of "CLimiterDetails". See corresponding hpp files for further details. + */ +template +void computeLimiters(ENUM_LIMITER LimiterKind, + CSolver* solver, + MPI_QUANTITIES kindMpiComm, + PERIODIC_QUANTITIES kindPeriodicComm1, + PERIODIC_QUANTITIES kindPeriodicComm2, + CGeometry& geometry, + CConfig& config, + size_t varBegin, + size_t varEnd, + const FieldType& field, + const GradientType& gradient, + FieldType& fieldMin, + FieldType& fieldMax, + FieldType& limiter) +{ +#define INSTANTIATE(KIND) \ +computeLimiters_impl(solver, kindMpiComm, \ + kindPeriodicComm1, kindPeriodicComm2, geometry, config, varBegin, \ + varEnd, field, gradient, fieldMin, fieldMax, limiter) + + switch (LimiterKind) { + case NO_LIMITER: + { + SU2_OMP_PARALLEL + { + SU2_OMP_FOR_STAT(512) + for(size_t iPoint = 0; iPoint < geometry.GetnPoint(); ++iPoint) + for(size_t iVar = varBegin; iVar < varEnd; ++iVar) + limiter(iPoint, iVar) = 1.0; + } + break; + } + case BARTH_JESPERSEN: + { + INSTANTIATE(BARTH_JESPERSEN); + break; + } + case VENKATAKRISHNAN: + { + INSTANTIATE(VENKATAKRISHNAN); + break; + } + case VENKATAKRISHNAN_WANG: + { + INSTANTIATE(VENKATAKRISHNAN_WANG); + break; + } + case WALL_DISTANCE: + { + INSTANTIATE(WALL_DISTANCE); + break; + } + case SHARP_EDGES: + { + INSTANTIATE(SHARP_EDGES); + break; + } + default: + { + SU2_MPI::Error("Unknown limiter type.", CURRENT_FUNCTION); + break; + } + } +#undef INSTANTIATE +} diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp new file mode 100644 index 000000000000..2c0b720c7034 --- /dev/null +++ b/SU2_CFD/include/limiters/computeLimiters_impl.hpp @@ -0,0 +1,260 @@ +/*! + * \file computeLimiters_impl.hpp + * \brief Generic and general computation of limiters. + * \note Common methods are derived by defining only small details. + * \author P. Gomes + * \version 7.0.0 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2019, 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 . + */ + + +/*! + * \brief General limiter computation for methods based on one limiter + * value per point (as opposed to one per edge) and per variable. + * \note This implementation can be used to derive most common methods + * by specializing the limiter functions (e.g. Venkatakrishnan) + * and the geometric modifications (e.g. sharp edges), this is done + * via specialization of "CLimiterDetails" (all its methods). Then, + * a call to this function with the specialized "LimiterKind" should + * be added to the body of "computeLimiters()". + * See also the notes in computeGradientsGreenGauss.hpp + * + * Arguments: + * \param[in] solver - Optional, solver associated with the field (used only for MPI). + * \param[in] kindMpiComm - Type of MPI communication required. + * \param[in] kindPeriodicComm1 - Type of periodic comm. to determine min/max. + * \param[in] kindPeriodicComm2 - Type of periodic comm. to adjust limiters. + * \param[in] geometry - Geometric grid properties. + * \param[in] config - Configuration of the problem. + * \param[in] varBegin - First variable index for which to compute limiters. + * \param[in] varEnd - End of computation range (nVar = end-begin). + * \param[in] field - Variable field. + * \param[in] gradient - Gradient of the field. + * \param[out] fieldMin - Minimum field values over direct neighbors of each point. + * \param[out] fieldMax - As above but maximum values. + * \param[out] limiter - Reconstruction limiter for the field. + * + * Template parameters: + * \param FieldType - Generic object with operator (iPoint,iVar) + * \param GradientType - Generic object with operator (iPoint,iVar,iDim) + * \param LimiterKind - Used to instantiate the right details class. + */ +template +void computeLimiters_impl(CSolver* solver, + MPI_QUANTITIES kindMpiComm, + PERIODIC_QUANTITIES kindPeriodicComm1, + PERIODIC_QUANTITIES kindPeriodicComm2, + CGeometry& geometry, + CConfig& config, + size_t varBegin, + size_t varEnd, + const FieldType& field, + const GradientType& gradient, + FieldType& fieldMin, + FieldType& fieldMax, + FieldType& limiter) +{ + constexpr size_t MAXNDIM = 3; + constexpr size_t MAXNVAR = 8; + + if (varEnd > MAXNVAR) + SU2_MPI::Error("Number of variables is too large, increase MAXNVAR.", CURRENT_FUNCTION); + + CLimiterDetails limiterDetails; + + size_t nPointDomain = geometry.GetnPointDomain(); + size_t nPoint = geometry.GetnPoint(); + size_t nDim = geometry.GetnDim(); + + /*--- If we do not have periodicity we can use a + * more efficient access pattern to memory. ---*/ + + bool periodic = (solver != nullptr) && + (kindPeriodicComm1 != PERIODIC_NONE) && + (config.GetnMarker_Periodic() > 0); + +#ifdef HAVE_OMP + constexpr size_t OMP_MAX_CHUNK = 512; + + size_t chunkSize = computeStaticChunkSize(nPointDomain, + omp_get_max_threads(), OMP_MAX_CHUNK); +#endif + +#ifdef CODI_REVERSE_TYPE + bool tapeActive = false; + + if (config.GetDiscrete_Adjoint() && config.GetFrozen_Limiter_Disc()) { + /*--- If limiters are frozen do not record the computation ---*/ + tapeActive = AD::globalTape.isActive(); + AD::StopRecording(); + } +#endif + + /*--- Start OpenMP parallel section. ---*/ + + SU2_OMP_PARALLEL + { + limiterDetails.preprocess(geometry, config, varBegin, varEnd, field); + + /*--- Courtesy barrier in case someone forgets preprocess is parallel. ---*/ + SU2_OMP_BARRIER + + /*--- Initialize all min/max field values if we have + * periodic comms. otherwise do it inside main loop. ---*/ + + if (periodic) + { + SU2_OMP_FOR_STAT(chunkSize) + for (size_t iPoint = 0; iPoint < nPoint; ++iPoint) + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + fieldMax(iPoint,iVar) = fieldMin(iPoint,iVar) = field(iPoint,iVar); + + SU2_OMP_MASTER + { + for (size_t iPeriodic = 1; iPeriodic <= config.GetnMarker_Periodic()/2; ++iPeriodic) + { + solver->InitiatePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm1); + solver->CompletePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm1); + } + } + SU2_OMP_BARRIER + } + + /*--- Compute limiter for each point. ---*/ + + SU2_OMP_FOR_DYN(chunkSize) + for (size_t iPoint = 0; iPoint < nPointDomain; ++iPoint) + { + auto node = geometry.node[iPoint]; + const su2double* coord_i = node->GetCoord(); + + AD::StartPreacc(); + AD::SetPreaccIn(coord_i, nDim); + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + { + AD::SetPreaccIn(field(iPoint,iVar)); + + if (periodic) { + /*--- Started outside loop, so counts as input. ---*/ + AD::SetPreaccIn(fieldMax(iPoint,iVar)); + AD::SetPreaccIn(fieldMin(iPoint,iVar)); + } + else { + /*--- Initialize min/max now for iPoint if not periodic. ---*/ + fieldMax(iPoint,iVar) = field(iPoint,iVar); + fieldMin(iPoint,iVar) = field(iPoint,iVar); + } + + for(size_t iDim = 0; iDim < nDim; ++iDim) + AD::SetPreaccIn(gradient(iPoint,iVar,iDim)); + } + + /*--- Initialize min/max projection out of iPoint. ---*/ + + su2double projMax[MAXNVAR], projMin[MAXNVAR]; + + for (size_t iVar = varBegin; iVar < varEnd; ++iVar) + projMax[iVar] = projMin[iVar] = 0.0; + + /*--- Compute max/min projection and values over direct neighbors. ---*/ + + for(size_t iNeigh = 0; iNeigh < node->GetnPoint(); ++iNeigh) + { + size_t jPoint = node->GetPoint(iNeigh); + + const su2double* coord_j = geometry.node[jPoint]->GetCoord(); + AD::SetPreaccIn(coord_j, nDim); + + /*--- Distance vector from iPoint to face (middle of the edge). ---*/ + + su2double dist_ij[MAXNDIM]; + + for(size_t iDim = 0; iDim < nDim; ++iDim) + dist_ij[iDim] = 0.5 * (coord_j[iDim] - coord_i[iDim]); + + /*--- Project each variable, update min/max. ---*/ + + for(size_t iVar = varBegin; iVar < varEnd; ++iVar) + { + su2double proj = 0.0; + + for(size_t iDim = 0; iDim < nDim; ++iDim) + proj += dist_ij[iDim] * gradient(iPoint,iVar,iDim); + + projMax[iVar] = max(projMax[iVar], proj); + projMin[iVar] = min(projMin[iVar], proj); + + AD::SetPreaccIn(field(jPoint,iVar)); + + fieldMax(iPoint,iVar) = max(fieldMax(iPoint,iVar), field(jPoint,iVar)); + fieldMin(iPoint,iVar) = min(fieldMin(iPoint,iVar), field(jPoint,iVar)); + } + } + + /*--- Compute the geometric factor. ---*/ + + su2double geoFactor = limiterDetails.geometricFactor(iPoint, geometry); + + /*--- Final limiter computation for each variable, get the min limiter + * out of the positive/negative projections and deltas. ---*/ + + for(size_t iVar = varBegin; iVar < varEnd; ++iVar) + { + su2double limMax = limiterDetails.limiterFunction(iVar, projMax[iVar], + fieldMax(iPoint,iVar) - field(iPoint,iVar)); + + su2double limMin = limiterDetails.limiterFunction(iVar, projMin[iVar], + fieldMin(iPoint,iVar) - field(iPoint,iVar)); + + limiter(iPoint,iVar) = geoFactor * min(limMax, limMin); + + AD::SetPreaccOut(limiter(iPoint,iVar)); + } + + AD::EndPreacc(); + } + + } // end SU2_OMP_PARALLEL + + /*--- If no solver was provided we do not communicate. ---*/ + + if (solver != nullptr) + { + /*--- Account for periodic effects, take the minimum limiter on each periodic pair. ---*/ + + for (size_t iPeriodic = 1; iPeriodic <= config.GetnMarker_Periodic()/2; ++iPeriodic) + { + solver->InitiatePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm2); + solver->CompletePeriodicComms(&geometry, &config, iPeriodic, kindPeriodicComm2); + } + + /*--- Obtain the limiters at halo points from the MPI ranks that own them. ---*/ + + solver->InitiateComms(&geometry, &config, kindMpiComm); + solver->CompleteComms(&geometry, &config, kindMpiComm); + } + +#ifdef CODI_REVERSE_TYPE + if (tapeActive) AD::StartRecording(); +#endif +} diff --git a/SU2_CFD/include/variables/CAdjEulerVariable.hpp b/SU2_CFD/include/variables/CAdjEulerVariable.hpp index e6292572450b..df96370ece95 100644 --- a/SU2_CFD/include/variables/CAdjEulerVariable.hpp +++ b/SU2_CFD/include/variables/CAdjEulerVariable.hpp @@ -167,4 +167,10 @@ class CAdjEulerVariable : public CVariable { */ inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } + /*! + * \brief Get the reconstruction gradient for variables at all points. + * \return Reference to reconstruction gradient. + */ + inline VectorOfMatrix& GetGradient_Reconstruction(void) final { return Gradient_Reconstruction; } + }; diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 8e2204cb2191..1842f57c5e73 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.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-2019, SU2 Contributors (cf. AUTHORS.md) @@ -94,11 +94,6 @@ class CEulerVariable : public CVariable { Solution_New(iPoint,iVar) += val_solution; } - /*! - * \brief Set to zero the gradient of the primitive variables. - */ - void SetGradient_PrimitiveZero() final; - /*! * \brief Add value to the gradient of the primitive variables. * \param[in] iVar - Index of the variable. @@ -109,16 +104,6 @@ class CEulerVariable : public CVariable { Gradient_Primitive(iPoint,iVar,iDim) += value; } - /*! - * \brief Subtract value to the gradient of the primitive variables. - * \param[in] iVar - Index of the variable. - * \param[in] iDim - Index of the dimension. - * \param[in] value - Value to subtract to the gradient of the primitive variables. - */ - inline void SubtractGradient_Primitive(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) final { - Gradient_Primitive(iPoint,iVar,iDim) -= value; - } - /*! * \brief Get the value of the primitive variables gradient. * \param[in] iVar - Index of the variable. @@ -129,6 +114,12 @@ class CEulerVariable : public CVariable { return Gradient_Primitive(iPoint,iVar,iDim); } + /*! + * \brief Get the primitive variables limiter. + * \return Primitive variables limiter for the entire domain. + */ + inline MatrixType& GetLimiter_Primitive(void) {return Limiter_Primitive; } + /*! * \brief Get the value of the primitive variables gradient. * \param[in] iVar - Index of the variable. @@ -155,6 +146,18 @@ class CEulerVariable : public CVariable { Limiter_Primitive(iPoint,iVar) = value; } + /*! + * \brief Get the primitive variable gradients for all points. + * \return Reference to primitive variable gradient. + */ + inline VectorOfMatrix& GetGradient_Primitive(void) { return Gradient_Primitive; } + + /*! + * \brief Get the reconstruction gradient for primitive variable at all points. + * \return Reference to variable reconstruction gradient. + */ + inline VectorOfMatrix& GetGradient_Reconstruction(void) final { return Gradient_Reconstruction; } + /*! * \brief Get the value of the primitive variables gradient. * \return Value of the primitive variables gradient. @@ -177,7 +180,7 @@ class CEulerVariable : public CVariable { inline su2double GetGradient_Reconstruction(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const final { return Gradient_Reconstruction(iPoint,iVar,iDim); } - + /*! * \brief Get the value of the reconstruction variables gradient at a node. * \param[in] iPoint - Index of the current node. @@ -188,14 +191,14 @@ class CEulerVariable : public CVariable { inline void SetGradient_Reconstruction(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) final { Gradient_Reconstruction(iPoint,iVar,iDim) = value; } - + /*! * \brief Get the array of the reconstruction variables gradient at a node. * \param[in] iPoint - Index of the current node. * \return Array of the reconstruction variables gradient at a node. */ inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } - + /*! * \brief A virtual member. */ @@ -244,6 +247,12 @@ class CEulerVariable : public CVariable { */ void SetSecondaryVar(unsigned long iPoint, CFluidModel *FluidModel); + /*! + * \brief Get the primitive variables for all points. + * \return Reference to primitives. + */ + inline const MatrixType& GetPrimitive(void) const { return Primitive; } + /*! * \brief Get the primitive variables. * \param[in] iVar - Index of the variable. diff --git a/SU2_CFD/include/variables/CHeatFVMVariable.hpp b/SU2_CFD/include/variables/CHeatFVMVariable.hpp index b96f95f97d1f..4c9bfe9e9626 100644 --- a/SU2_CFD/include/variables/CHeatFVMVariable.hpp +++ b/SU2_CFD/include/variables/CHeatFVMVariable.hpp @@ -87,4 +87,10 @@ class CHeatFVMVariable final : public CVariable { */ inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } + /*! + * \brief Get the reconstruction gradient for primitive variable at all points. + * \return Reference to variable reconstruction gradient. + */ + inline VectorOfMatrix& GetGradient_Reconstruction(void) final { return Gradient_Reconstruction; } + }; diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 742089b08394..712daf6ced76 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -65,9 +65,16 @@ class CIncEulerVariable : public CVariable { virtual ~CIncEulerVariable() = default; /*! - * \brief Set to zero the gradient of the primitive variables. + * \brief Get the primitive variable gradients for all points. + * \return Reference to primitive variable gradient. */ - void SetGradient_PrimitiveZero() final; + inline VectorOfMatrix& GetGradient_Primitive(void) { return Gradient_Primitive; } + + /*! + * \brief Get the reconstruction gradient for primitive variable at all points. + * \return Reference to variable reconstruction gradient. + */ + inline VectorOfMatrix& GetGradient_Reconstruction(void) final { return Gradient_Reconstruction; } /*! * \brief Add value to the gradient of the primitive variables. @@ -80,17 +87,6 @@ class CIncEulerVariable : public CVariable { Gradient_Primitive(iPoint,iVar,iDim) += value; } - /*! - * \brief Subtract value to the gradient of the primitive variables. - * \param[in] iPoint - Point index. - * \param[in] iVar - Index of the variable. - * \param[in] iDim - Index of the dimension. - * \param[in] value - Value to subtract to the gradient of the primitive variables. - */ - inline void SubtractGradient_Primitive(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) final { - Gradient_Primitive(iPoint,iVar,iDim) -= value; - } - /*! * \brief Get the value of the primitive variables gradient. * \param[in] iPoint - Point index. @@ -102,6 +98,12 @@ class CIncEulerVariable : public CVariable { return Gradient_Primitive(iPoint,iVar,iDim); } + /*! + * \brief Get the primitive variables limiter. + * \return Primitive variables limiter for the entire domain. + */ + inline MatrixType& GetLimiter_Primitive(void) {return Limiter_Primitive; } + /*! * \brief Get the value of the primitive variables gradient. * \param[in] iPoint - Point index. @@ -182,6 +184,12 @@ class CIncEulerVariable : public CVariable { */ inline void SetPressure(unsigned long iPoint) final { Primitive(iPoint,0) = Solution(iPoint,0); } + /*! + * \brief Get the primitive variables for all points. + * \return Reference to primitives. + */ + inline const MatrixType& GetPrimitive(void) const { return Primitive; } + /*! * \brief Get the primitive variables. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index 48a232997938..d69659bc9070 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -100,6 +100,12 @@ class CTurbVariable : public CVariable { * \return Array of the reconstruction variables gradient at a node. */ inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } - + + /*! + * \brief Get the reconstruction gradient for primitive variable at all points. + * \return Reference to variable reconstruction gradient. + */ + inline VectorOfMatrix& GetGradient_Reconstruction(void) final { return Gradient_Reconstruction; } + }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index cdb3301c55c6..afa55eb3a38d 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.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-2019, SU2 Contributors (cf. AUTHORS.md) @@ -82,11 +82,13 @@ class CVariable { MatrixType External; /*!< \brief External (outer) contribution in discrete adjoint multizone problems. */ su2vector Non_Physical; /*!< \brief Non-physical points in the solution (force first order). */ - su2vector Non_Physical_Counter; /*!< \brief Number of consecutive iterations that a point has been treated first-order. After a specified number of successful reconstructions, the point can be returned to second-order. */ - + su2vector + Non_Physical_Counter; /*!< \brief Number of consecutive iterations that a point has been treated first-order. + After a specified number of successful reconstructions, the point can be returned to second-order. */ + VectorType UnderRelaxation; /*!< \brief Value of the under-relxation parameter local to the control volume. */ VectorType LocalCFL; /*!< \brief Value of the CFL number local to the control volume. */ - + MatrixType Solution_time_n; /*!< \brief Solution of the problem at time n for dual-time stepping technique. */ MatrixType Solution_time_n1; /*!< \brief Solution of the problem at time n-1 for dual-time stepping technique. */ VectorType Delta_Time; /*!< \brief Time step. */ @@ -458,6 +460,12 @@ class CVariable { Solution(iPoint,iVar) = min(max(val_new, lowerlimit), upperlimit); } + /*! + * \brief Get the entire solution of the problem. + * \return Reference to the solution matrix. + */ + inline const MatrixType& GetSolution(void) { return Solution; } + /*! * \brief Get the solution of the problem. * \param[in] iPoint - Point index. @@ -554,21 +562,21 @@ class CVariable { * \param[in] val_under_relaxation - the input value of the under-relaxation parameter for this CV. */ inline void SetUnderRelaxation(unsigned long iPoint, su2double val_under_relaxation) { UnderRelaxation(iPoint) = val_under_relaxation; } - + /*! * \brief Get the value of the under-relaxation parameter for the current control volume (CV). * \param[in] iPoint - Point index. * \return Value of the under-relaxation parameter for this CV. */ inline su2double GetUnderRelaxation(unsigned long iPoint) const { return UnderRelaxation(iPoint); } - + /*! * \brief Set the value of the local CFL number for the current control volume (CV). * \param[in] iPoint - Point index. * \param[in] val_cfl - the input value of the local CFL number for this CV. */ inline void SetLocalCFL(unsigned long iPoint, su2double val_cfl) { LocalCFL(iPoint) = val_cfl; } - + /*! * \brief Get the value of the local CFL number for the current control volume (CV). * \param[in] iPoint - Point index. @@ -591,9 +599,12 @@ class CVariable { inline su2double GetAuxVar(unsigned long iPoint) const { return AuxVar(iPoint); } /*! - * \brief Set the auxiliary variable gradient to zero value. + * \brief Get the auxiliary variable. + * \return 2D view of the auxiliary variable. */ - void SetAuxVarGradientZero(); + inline C2DDummyLastView GetAuxVar(void) const { + return C2DDummyLastView(AuxVar); + } /*! * \brief Set the value of the auxiliary variable gradient. @@ -612,19 +623,19 @@ class CVariable { inline void AddAuxVarGradient(unsigned long iPoint, unsigned long iDim, su2double val_value) { Grad_AuxVar(iPoint,iDim) += val_value;} /*! - * \brief Subtract a value to the auxiliary variable gradient. + * \brief Get the gradient of the auxiliary variable. * \param[in] iPoint - Point index. - * \param[in] iDim - Index of the dimension. - * \param[in] val_value - Value of the gradient to be subtracted for the index iDim. + * \return Value of the gradient of the auxiliary variable. */ - inline void SubtractAuxVarGradient(unsigned long iPoint, unsigned long iDim, su2double val_value) { Grad_AuxVar(iPoint,iDim) -= val_value; } + inline su2double *GetAuxVarGradient(unsigned long iPoint) { return Grad_AuxVar[iPoint]; } /*! * \brief Get the gradient of the auxiliary variable. - * \param[in] iPoint - Point index. - * \return Value of the gradient of the auxiliary variable. + * \return 3D view of the gradient of the auxiliary variable. */ - inline su2double *GetAuxVarGradient(unsigned long iPoint) { return Grad_AuxVar[iPoint]; } + inline C3DDummyMiddleView GetAuxVarGradient() { + return C3DDummyMiddleView(Grad_AuxVar); + } /*! * \brief Get the gradient of the auxiliary variable. @@ -719,11 +730,6 @@ class CVariable { */ inline void SetGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) { Gradient(iPoint,iVar,iDim) = value; } - /*! - * \brief Set to zero the gradient of the solution. - */ - void SetGradientZero(); - /*! * \brief Add value to the solution gradient. * \param[in] iPoint - Point index. @@ -734,13 +740,10 @@ class CVariable { inline void AddGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) { Gradient(iPoint,iVar,iDim) += value; } /*! - * \brief Subtract value to the solution gradient. - * \param[in] iPoint - Point index. - * \param[in] iVar - Index of the variable. - * \param[in] iDim - Index of the dimension. - * \param[in] value - Value to subtract to the solution gradient. + * \brief Get the gradient of the entire solution. + * \return Reference to gradient. */ - inline void SubtractGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) { Gradient(iPoint,iVar,iDim) -= value; } + inline VectorOfMatrix& GetGradient(void) { return Gradient; } /*! * \brief Get the value of the solution gradient. @@ -758,20 +761,6 @@ class CVariable { */ inline su2double GetGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return Gradient(iPoint,iVar,iDim); } - /*! - * \brief Set the value of an entry in the Rmatrix for least squares gradient calculations. - * \param[in] iPoint - Point index. - * \param[in] iDim - Index of the dimension. - * \param[in] jDim - Index of the dimension. - * \param[in] value - Value of the Rmatrix entry. - */ - inline void SetRmatrix(unsigned long iPoint, unsigned long iDim, unsigned long jDim, su2double value) { Rmatrix(iPoint,iDim,jDim) = value; } - - /*! - * \brief Set to zero the Rmatrix for least squares gradient calculations. - */ - void SetRmatrixZero(); - /*! * \brief Add value to the Rmatrix for least squares gradient calculations. * \param[in] iPoint - Point index. @@ -797,6 +786,12 @@ class CVariable { */ inline su2double **GetRmatrix(unsigned long iPoint) { return Rmatrix[iPoint]; } + /*! + * \brief Get the value Rmatrix for the entire domain. + * \return Reference to the Rmatrix. + */ + inline VectorOfMatrix& GetRmatrix(void) { return Rmatrix; } + /*! * \brief Set the value of the limiter. * \param[in] iPoint - Point index. @@ -838,6 +833,12 @@ class CVariable { */ inline void SetSolution_Min(unsigned long iPoint, unsigned long iVar, su2double solution) { Solution_Min(iPoint,iVar) = solution; } + /*! + * \brief Get the slope limiter. + * \return Reference to the limiters vector. + */ + inline MatrixType& GetLimiter(void) { return Limiter; } + /*! * \brief Get the value of the slope limiter. * \param[in] iPoint - Point index. @@ -861,6 +862,12 @@ class CVariable { */ inline su2double GetSolution_Max(unsigned long iPoint, unsigned long iVar) const { return Solution_Max(iPoint,iVar); } + /*! + * \brief Get the min solution. + * \return Value of the min solution for the domain. + */ + inline MatrixType& GetSolution_Max(void) { return Solution_Max; } + /*! * \brief Set the value of the preconditioner Beta. * \param[in] val_Beta - Value of the low Mach preconditioner variable Beta @@ -870,6 +877,12 @@ class CVariable { */ inline su2double GetSolution_Min(unsigned long iPoint, unsigned long iVar) const { return Solution_Min(iPoint,iVar); } + /*! + * \brief Get the min solution. + * \return Value of the min solution for the domain. + */ + inline MatrixType& GetSolution_Min(void) { return Solution_Min; } + /*! * \brief Get the value of the wind gust * \param[in] iPoint - Point index. @@ -1897,11 +1910,6 @@ class CVariable { */ inline virtual void SetVelSolutionDVector(unsigned long iPoint) {} - /*! - * \brief A virtual member. - */ - virtual void SetGradient_PrimitiveZero() {} - /*! * \brief A virtual member. * \param[in] iVar - Index of the variable. @@ -1910,14 +1918,6 @@ class CVariable { */ inline virtual void AddGradient_Primitive(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double val_value) {} - /*! - * \brief A virtual member. - * \param[in] iVar - Index of the variable. - * \param[in] iDim - Index of the dimension. - * \param[in] val_value - Value to subtract to the gradient of the primitive variables. - */ - inline virtual void SubtractGradient_Primitive(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double val_value) {} - /*! * \brief A virtual member. * \param[in] iVar - Index of the variable. @@ -1967,7 +1967,7 @@ class CVariable { * \return Value of the primitive variables gradient. */ inline virtual su2double GetGradient_Reconstruction(unsigned long iPoint, unsigned long val_var, unsigned long val_dim) const { return 0.0; } - + /*! * \brief Set the value of the primitive gradient for MUSCL reconstruction. * \param[in] val_var - Index of the variable. @@ -1975,13 +1975,19 @@ class CVariable { * \param[in] val_value - Value of the gradient. */ inline virtual void SetGradient_Reconstruction(unsigned long iPoint, unsigned long val_var, unsigned long val_dim, su2double val_value) {} - + /*! * \brief Get the value of the primitive gradient for MUSCL reconstruction. * \return Value of the primitive gradient for MUSCL reconstruction. */ inline virtual su2double **GetGradient_Reconstruction(unsigned long iPoint) { return nullptr; } - + + /*! + * \brief Get the reconstruction gradient for primitive variable at all points. + * \return Reference to variable reconstruction gradient. + */ + inline virtual VectorOfMatrix& GetGradient_Reconstruction(void) { return Gradient; } + /*! * \brief Set the blending function for the blending of k-w and k-eps. * \param[in] val_viscosity - Value of the vicosity. diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am index be09c7ecd791..11ac6b8cc98d 100644 --- a/SU2_CFD/obj/Makefile.am +++ b/SU2_CFD/obj/Makefile.am @@ -160,7 +160,8 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/variables/CAdjEulerVariable.cpp \ ../src/variables/CDiscAdjVariable.cpp \ ../src/variables/CIncNSVariable.cpp \ - ../src/variables/CEulerVariable.cpp + ../src/variables/CEulerVariable.cpp \ + ../src/limiters/computeLimiters.cpp su2_cfd_sources = \ ../include/SU2_CFD.hpp \ diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index fa341db4d593..455a95b5c83b 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -4178,6 +4178,19 @@ bool CFluidDriver::Monitor(unsigned long ExtIter) { void CFluidDriver::Output(unsigned long InnerIter) { + for (iZone = 0; iZone < nZone; iZone++) { + const auto inst = config_container[iZone]->GetiInst(); + + for (iInst = 0; iInst < nInst[iZone]; ++iInst) { + config_container[iZone]->SetiInst(iInst); + output_container[iZone]->SetResult_Files(geometry_container[iZone][iInst][MESH_0], + config_container[iZone], + solver_container[iZone][iInst][MESH_0], + InnerIter, StopCalc); + } + config_container[iZone]->SetiInst(inst); + } + } diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index 9ca16fd27b36..82360d43437d 100644 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -30,6 +30,10 @@ #include "../../Common/include/toolboxes/printing_toolbox.hpp" #include "../include/variables/CEulerVariable.hpp" #include "../include/variables/CNSVariable.hpp" +#include "../include/gradients/computeGradientsGreenGauss.hpp" +#include "../include/gradients/computeGradientsLeastSquares.hpp" +#include "../include/limiters/computeLimiters.hpp" + CEulerSolver::CEulerSolver(void) : CSolver() { @@ -5268,620 +5272,44 @@ void CEulerSolver::ComputeUnderRelaxationFactor(CSolver **solver_container, CCon } void CEulerSolver::SetPrimitive_Gradient_GG(CGeometry *geometry, CConfig *config, bool reconstruction) { - unsigned long iPoint, jPoint, iEdge, iVertex; - unsigned short iDim, iVar, iMarker; - su2double *PrimVar_Vertex, *PrimVar_i, *PrimVar_j, PrimVar_Average, - Partial_Gradient, Partial_Res, *Normal, Vol; - - /*--- Gradient primitive variables compressible (temp, vx, vy, vz, P, rho) ---*/ - - PrimVar_Vertex = new su2double [nPrimVarGrad]; - PrimVar_i = new su2double [nPrimVarGrad]; - PrimVar_j = new su2double [nPrimVarGrad]; - - /*--- Set Gradient_Primitive to zero ---*/ - - nodes->SetGradient_PrimitiveZero(); - - /*--- Loop interior edges ---*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - PrimVar_i[iVar] = nodes->GetPrimitive(iPoint,iVar); - PrimVar_j[iVar] = nodes->GetPrimitive(jPoint,iVar); - } - - Normal = geometry->edge[iEdge]->GetNormal(); - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - PrimVar_Average = 0.5 * ( PrimVar_i[iVar] + PrimVar_j[iVar] ); - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Res = PrimVar_Average*Normal[iDim]; - if (geometry->node[iPoint]->GetDomain()) - nodes->AddGradient_Primitive(iPoint,iVar, iDim, Partial_Res); - if (geometry->node[jPoint]->GetDomain()) - nodes->SubtractGradient_Primitive(jPoint,iVar, iDim, Partial_Res); - } - } - } - - /*--- Loop boundary edges ---*/ - - for (iMarker = 0; 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++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - if (geometry->node[iPoint]->GetDomain()) { - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - PrimVar_Vertex[iVar] = nodes->GetPrimitive(iPoint,iVar); - - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Res = PrimVar_Vertex[iVar]*Normal[iDim]; - nodes->SubtractGradient_Primitive(iPoint,iVar, iDim, Partial_Res); - } - } - } - } - } - /*--- Correct the sensor values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_GG); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_GG); - } - - /*--- Update gradient value ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Get the volume, which may include periodic components. ---*/ - - Vol = (geometry->node[iPoint]->GetVolume() + - geometry->node[iPoint]->GetPeriodicVolume()); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Gradient = nodes->GetGradient_Primitive(iPoint, iVar, iDim)/Vol; - if (reconstruction) - nodes->SetGradient_Reconstruction(iPoint, iVar, iDim, Partial_Gradient); - else - nodes->SetGradient_Primitive(iPoint, iVar, iDim, Partial_Gradient); - } - } - - } + const auto& primitives = nodes->GetPrimitive(); + auto& gradient = reconstruction? nodes->GetGradient_Reconstruction() : nodes->GetGradient_Primitive(); - delete [] PrimVar_Vertex; - delete [] PrimVar_i; - delete [] PrimVar_j; - - /*--- Communicate the gradient values via MPI. ---*/ - - InitiateComms(geometry, config, PRIMITIVE_GRADIENT); - CompleteComms(geometry, config, PRIMITIVE_GRADIENT); - + computeGradientsGreenGauss(this, PRIMITIVE_GRADIENT, PERIODIC_PRIM_GG, *geometry, + *config, primitives, 0, nPrimVarGrad, gradient); } void CEulerSolver::SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config, bool reconstruction) { - - unsigned short iVar, iDim, jDim, iNeigh; - unsigned long iPoint, jPoint; - su2double *PrimVar_i, *PrimVar_j, *Coord_i, *Coord_j; - su2double r11, r12, r13, r22, r23, r23_a, r23_b, r33, weight; - su2double z11, z12, z13, z22, z23, z33, detR2; - bool singular; - - /*--- Set a flag for unweighted or weighted least-squares. ---*/ - - bool weighted = true; - if (reconstruction) { - if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) - weighted = false; - } else if (config->GetKind_Gradient_Method() == LEAST_SQUARES) { - weighted = false; - } - /*--- Clear Rmatrix and Primitive gradient ---*/ - nodes->SetRmatrixZero(); - nodes->SetGradient_PrimitiveZero(); - - /*--- Loop over points of the grid ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Set the value of the singular ---*/ - singular = false; - - /*--- Get coordinates ---*/ - - Coord_i = geometry->node[iPoint]->GetCoord(); - - /*--- Get primitives from CVariable ---*/ - - PrimVar_i = nodes->GetPrimitive(iPoint); - - /*--- Inizialization of variables ---*/ - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - for (iDim = 0; iDim < nDim; iDim++) - Cvector[iVar][iDim] = 0.0; - - /*--- Clear Rmatrix, which could eventually be computed once - and stored for static meshes, as well as the prim gradient. ---*/ - - AD::StartPreacc(); - AD::SetPreaccIn(PrimVar_i, nPrimVarGrad); - AD::SetPreaccIn(Coord_i, nDim); - - for (iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); iNeigh++) { - jPoint = geometry->node[iPoint]->GetPoint(iNeigh); - Coord_j = geometry->node[jPoint]->GetCoord(); + /*--- Set a flag for unweighted or weighted least-squares. ---*/ + bool weighted; - PrimVar_j = nodes->GetPrimitive(jPoint); + if (reconstruction) + weighted = (config->GetKind_Gradient_Method_Recon() == WEIGHTED_LEAST_SQUARES); + else + weighted = (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES); - AD::SetPreaccIn(Coord_j, nDim); - AD::SetPreaccIn(PrimVar_j, nPrimVarGrad); - - if (weighted) { - weight = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - weight += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - } else { - weight = 1.0; - } - - /*--- Sumations for entries of upper triangular matrix R ---*/ - - if (weight != 0.0) { + const auto& primitives = nodes->GetPrimitive(); + auto& rmatrix = nodes->GetRmatrix(); + auto& gradient = reconstruction? nodes->GetGradient_Reconstruction() : nodes->GetGradient_Primitive(); + PERIODIC_QUANTITIES kindPeriodicComm = weighted? PERIODIC_PRIM_LS : PERIODIC_PRIM_ULS; - nodes->AddRmatrix(iPoint,0, 0, (Coord_j[0]-Coord_i[0])*(Coord_j[0]-Coord_i[0])/weight); - nodes->AddRmatrix(iPoint,0, 1, (Coord_j[0]-Coord_i[0])*(Coord_j[1]-Coord_i[1])/weight); - nodes->AddRmatrix(iPoint,1, 1, (Coord_j[1]-Coord_i[1])*(Coord_j[1]-Coord_i[1])/weight); - - if (nDim == 3) { - nodes->AddRmatrix(iPoint,0, 2, (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight); - nodes->AddRmatrix(iPoint,1, 2, (Coord_j[1]-Coord_i[1])*(Coord_j[2]-Coord_i[2])/weight); - nodes->AddRmatrix(iPoint,2, 1, (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight); - nodes->AddRmatrix(iPoint,2, 2, (Coord_j[2]-Coord_i[2])*(Coord_j[2]-Coord_i[2])/weight); - } - - /*--- Entries of c:= transpose(A)*b ---*/ - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - nodes->AddGradient_Primitive(iPoint,iVar,iDim, (Coord_j[iDim]-Coord_i[iDim])*(PrimVar_j[iVar]-PrimVar_i[iVar])/weight); - } - } - } - } - AD::SetPreaccOut(nodes->GetRmatrix(iPoint), nDim, nDim); - AD::SetPreaccOut(nodes->GetGradient_Primitive(iPoint), nPrimVarGrad, nDim); - AD::EndPreacc(); - } - - /*--- Correct the gradient values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - if (weighted) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_LS); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_LS); - } else { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_ULS); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_ULS); - } - } - - /*--- Second loop over points of the grid to compute final gradient ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Set the value of the singular ---*/ - - singular = false; - - /*--- Entries of upper triangular matrix R ---*/ - - r11 = 0.0; r12 = 0.0; r13 = 0.0; r22 = 0.0; - r23 = 0.0; r23_a = 0.0; r23_b = 0.0; r33 = 0.0; - - r11 = nodes->GetRmatrix(iPoint,0,0); - r12 = nodes->GetRmatrix(iPoint,0,1); - r22 = nodes->GetRmatrix(iPoint,1,1); - - AD::StartPreacc(); - AD::SetPreaccIn(r11); - AD::SetPreaccIn(r12); - AD::SetPreaccIn(r22); - - if (r11 >= 0.0) r11 = sqrt(r11); else r11 = 0.0; - if (r11 != 0.0) r12 = r12/r11; else r12 = 0.0; - if (r22-r12*r12 >= 0.0) r22 = sqrt(r22-r12*r12); else r22 = 0.0; - - if (nDim == 3) { - r13 = nodes->GetRmatrix(iPoint,0,2); - r23_a = nodes->GetRmatrix(iPoint,1,2); - r23_b = nodes->GetRmatrix(iPoint,2,1); - r33 = nodes->GetRmatrix(iPoint,2,2); - - AD::SetPreaccIn(r13); - AD::SetPreaccIn(r23_a); - AD::SetPreaccIn(r23_b); - AD::SetPreaccIn(r33); - - if (r11 != 0.0) r13 = r13/r11; else r13 = 0.0; - if ((r22 != 0.0) && (r11*r22 != 0.0)) r23 = r23_a/r22 - r23_b*r12/(r11*r22); else r23 = 0.0; - if (r33-r23*r23-r13*r13 >= 0.0) r33 = sqrt(r33-r23*r23-r13*r13); else r33 = 0.0; - } - - /*--- Compute determinant ---*/ - - if (nDim == 2) detR2 = (r11*r22)*(r11*r22); - else detR2 = (r11*r22*r33)*(r11*r22*r33); - - /*--- Detect singular matrices ---*/ - - if (abs(detR2) <= EPS) { detR2 = 1.0; singular = true; } - - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - - if (singular) { - for (iDim = 0; iDim < nDim; iDim++) - for (jDim = 0; jDim < nDim; jDim++) - Smatrix[iDim][jDim] = 0.0; - } - else { - if (nDim == 2) { - Smatrix[0][0] = (r12*r12+r22*r22)/detR2; - Smatrix[0][1] = -r11*r12/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = r11*r11/detR2; - } - else { - z11 = r22*r33; z12 = -r12*r33; z13 = r12*r23-r13*r22; - z22 = r11*r33; z23 = -r11*r23; z33 = r11*r22; - Smatrix[0][0] = (z11*z11+z12*z12+z13*z13)/detR2; - Smatrix[0][1] = (z12*z22+z13*z23)/detR2; - Smatrix[0][2] = (z13*z33)/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = (z22*z22+z23*z23)/detR2; - Smatrix[1][2] = (z23*z33)/detR2; - Smatrix[2][0] = Smatrix[0][2]; - Smatrix[2][1] = Smatrix[1][2]; - Smatrix[2][2] = (z33*z33)/detR2; - } - } - - AD::SetPreaccOut(Smatrix, nDim, nDim); - AD::EndPreacc(); - - /*--- Computation of the gradient: S*c ---*/ - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - Cvector[iVar][iDim] = 0.0; - for (jDim = 0; jDim < nDim; jDim++) { - Cvector[iVar][iDim] += Smatrix[iDim][jDim]*nodes->GetGradient_Primitive(iPoint,iVar, jDim); - } - } - } - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - if (reconstruction) - nodes->SetGradient_Reconstruction(iPoint, iVar, iDim, Cvector[iVar][iDim]); - else - nodes->SetGradient_Primitive(iPoint, iVar, iDim, Cvector[iVar][iDim]); - } - } - - } - - /*--- Communicate the gradient values via MPI. ---*/ - - InitiateComms(geometry, config, PRIMITIVE_GRADIENT); - CompleteComms(geometry, config, PRIMITIVE_GRADIENT); - + computeGradientsLeastSquares(this, PRIMITIVE_GRADIENT, kindPeriodicComm, *geometry, *config, + weighted, primitives, 0, nPrimVarGrad, gradient, rmatrix); } void CEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) { - - unsigned long iEdge, iPoint, jPoint; - unsigned short iVar, iDim; - su2double **Gradient_i, **Gradient_j, *Coord_i, *Coord_j, - *Primitive, *Primitive_i, *Primitive_j, - *LocalMinPrimitive = NULL, *LocalMaxPrimitive = NULL, - *GlobalMinPrimitive = NULL, *GlobalMaxPrimitive = NULL, - dave, LimK, eps2, eps1, dm, dp, du, y, limiter; - -#ifdef CODI_REVERSE_TYPE - bool TapeActive = false; - - if (config->GetDiscrete_Adjoint() && config->GetFrozen_Limiter_Disc()) { - /*--- If limiters are frozen do not record the computation ---*/ - TapeActive = AD::globalTape.isActive(); - AD::StopRecording(); - } -#endif - - dave = config->GetRefElemLength(); - LimK = config->GetVenkat_LimiterCoeff(); - - if (config->GetKind_SlopeLimit_Flow() == NO_LIMITER) { - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - nodes->SetLimiter_Primitive(iPoint,iVar, 1.0); - } - } - - } - - else { - - /*--- Initialize solution max and solution min and the limiter in the entire domain --*/ - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - nodes->SetSolution_Max(iPoint,iVar, -EPS); - nodes->SetSolution_Min(iPoint,iVar, EPS); - nodes->SetLimiter_Primitive(iPoint,iVar, 2.0); - } - } - - /*--- Establish bounds for Spekreijse monotonicity by finding max & min values of neighbor variables --*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - - /*--- Get the primitive variables ---*/ - - Primitive_i = nodes->GetPrimitive(iPoint); - Primitive_j = nodes->GetPrimitive(jPoint); - - /*--- Compute the maximum, and minimum values for nodes i & j ---*/ - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - du = (Primitive_j[iVar] - Primitive_i[iVar]); - nodes->SetSolution_Min(iPoint, iVar, min(nodes->GetSolution_Min(iPoint, iVar), du)); - nodes->SetSolution_Max(iPoint, iVar, max(nodes->GetSolution_Max(iPoint, iVar), du)); - nodes->SetSolution_Min(jPoint, iVar, min(nodes->GetSolution_Min(jPoint, iVar), -du)); - nodes->SetSolution_Max(jPoint, iVar, max(nodes->GetSolution_Max(jPoint, iVar), -du)); - } - - } - - /*--- Correct the limiter values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_PRIM_1); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_PRIM_1); - } - - } - - /*--- Barth-Jespersen limiter with Venkatakrishnan modification ---*/ - - if (config->GetKind_SlopeLimit_Flow() == BARTH_JESPERSEN) { - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Gradient_i = nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = nodes->GetGradient_Reconstruction(jPoint); - Coord_i = geometry->node[iPoint]->GetCoord(); - Coord_j = geometry->node[jPoint]->GetCoord(); - - AD::StartPreacc(); - AD::SetPreaccIn(Gradient_i, nPrimVarGrad, nDim); - AD::SetPreaccIn(Gradient_j, nPrimVarGrad, nDim); - AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - - AD::SetPreaccIn(nodes->GetSolution_Max(iPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Min(iPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Max(jPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Min(jPoint,iVar)); - - /*--- Calculate the interface left gradient, delta- (dm) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_j[iDim]-Coord_i[iDim])*Gradient_i[iVar][iDim]; - - if (dm == 0.0) { limiter = 2.0; } - else { - if ( dm > 0.0 ) dp = nodes->GetSolution_Max(iPoint,iVar); - else dp = nodes->GetSolution_Min(iPoint,iVar); - limiter = dp/dm; - } - - if (limiter < nodes->GetLimiter_Primitive(iPoint,iVar)) { - nodes->SetLimiter_Primitive(iPoint,iVar, limiter); - AD::SetPreaccOut(nodes->GetLimiter_Primitive(iPoint)[iVar]); - } - - /*--- Calculate the interface right gradient, delta+ (dp) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_i[iDim]-Coord_j[iDim])*Gradient_j[iVar][iDim]; - - if (dm == 0.0) { limiter = 2.0; } - else { - if ( dm > 0.0 ) dp = nodes->GetSolution_Max(jPoint,iVar); - else dp = nodes->GetSolution_Min(jPoint,iVar); - limiter = dp/dm; - } - - if (limiter < nodes->GetLimiter_Primitive(jPoint,iVar)) { - nodes->SetLimiter_Primitive(jPoint,iVar, limiter); - AD::SetPreaccOut(nodes->GetLimiter_Primitive(jPoint)[iVar]); - } - - } - - AD::EndPreacc(); - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - y = nodes->GetLimiter_Primitive(iPoint,iVar); - limiter = (y*y + 2.0*y) / (y*y + y + 2.0); - nodes->SetLimiter_Primitive(iPoint,iVar, limiter); - } - } - - } - - /*--- Venkatakrishnan limiter ---*/ - - if ((config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN) || - (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG)) { - - if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { + auto kindLimiter = static_cast(config->GetKind_SlopeLimit_Flow()); + const auto& primitives = nodes->GetPrimitive(); + const auto& gradient = nodes->GetGradient_Reconstruction(); + auto& primMin = nodes->GetSolution_Min(); + auto& primMax = nodes->GetSolution_Max(); + auto& limiter = nodes->GetLimiter_Primitive(); - /*--- Allocate memory for the max and min primitive value --*/ - - LocalMinPrimitive = new su2double [nPrimVarGrad]; GlobalMinPrimitive = new su2double [nPrimVarGrad]; - LocalMaxPrimitive = new su2double [nPrimVarGrad]; GlobalMaxPrimitive = new su2double [nPrimVarGrad]; - - /*--- Compute the max value and min value of the solution ---*/ - - Primitive = nodes->GetPrimitive(0); - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - LocalMinPrimitive[iVar] = Primitive[iVar]; - LocalMaxPrimitive[iVar] = Primitive[iVar]; - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - - /*--- Get the primitive variables ---*/ - - Primitive = nodes->GetPrimitive(iPoint); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - LocalMinPrimitive[iVar] = min (LocalMinPrimitive[iVar], Primitive[iVar]); - LocalMaxPrimitive[iVar] = max (LocalMaxPrimitive[iVar], Primitive[iVar]); - } - - } - -#ifdef HAVE_MPI - SU2_MPI::Allreduce(LocalMinPrimitive, GlobalMinPrimitive, nPrimVarGrad, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - SU2_MPI::Allreduce(LocalMaxPrimitive, GlobalMaxPrimitive, nPrimVarGrad, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); -#else - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - GlobalMinPrimitive[iVar] = LocalMinPrimitive[iVar]; - GlobalMaxPrimitive[iVar] = LocalMaxPrimitive[iVar]; - } -#endif - } - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Gradient_i = nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = nodes->GetGradient_Reconstruction(jPoint); - Coord_i = geometry->node[iPoint]->GetCoord(); - Coord_j = geometry->node[jPoint]->GetCoord(); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - - AD::StartPreacc(); - AD::SetPreaccIn(Gradient_i[iVar], nDim); - AD::SetPreaccIn(Gradient_j[iVar], nDim); - AD::SetPreaccIn(Coord_i, nDim); - AD::SetPreaccIn(Coord_j, nDim); - AD::SetPreaccIn(nodes->GetSolution_Max(iPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Min(iPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Max(jPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Min(jPoint,iVar)); - - if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { - AD::SetPreaccIn(GlobalMaxPrimitive[iVar]); - AD::SetPreaccIn(GlobalMinPrimitive[iVar]); - eps1 = LimK * (GlobalMaxPrimitive[iVar] - GlobalMinPrimitive[iVar]); - eps2 = eps1*eps1; - } - else { - eps1 = LimK*dave; - eps2 = eps1*eps1*eps1; - } - - /*--- Calculate the interface left gradient, delta- (dm) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_j[iDim]-Coord_i[iDim])*Gradient_i[iVar][iDim]; - - /*--- Calculate the interface right gradient, delta+ (dp) ---*/ - - if ( dm > 0.0 ) dp = nodes->GetSolution_Max(iPoint,iVar); - else dp = nodes->GetSolution_Min(iPoint,iVar); - - limiter = ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < nodes->GetLimiter_Primitive(iPoint,iVar)) { - nodes->SetLimiter_Primitive(iPoint,iVar, limiter); - AD::SetPreaccOut(nodes->GetLimiter_Primitive(iPoint)[iVar]); - } - - /*-- Repeat for point j on the edge ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_i[iDim]-Coord_j[iDim])*Gradient_j[iVar][iDim]; - - if ( dm > 0.0 ) dp = nodes->GetSolution_Max(jPoint,iVar); - else dp = nodes->GetSolution_Min(jPoint,iVar); - - limiter = ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < nodes->GetLimiter_Primitive(jPoint,iVar)) { - nodes->SetLimiter_Primitive(jPoint,iVar, limiter); - AD::SetPreaccOut(nodes->GetLimiter_Primitive(jPoint)[iVar]); - } - - AD::EndPreacc(); - } - } - - if (LocalMinPrimitive != NULL) delete [] LocalMinPrimitive; - if (LocalMaxPrimitive != NULL) delete [] LocalMaxPrimitive; - if (GlobalMinPrimitive != NULL) delete [] GlobalMinPrimitive; - if (GlobalMaxPrimitive != NULL) delete [] GlobalMaxPrimitive; - - } - - /*--- Correct the limiter values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_PRIM_2); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_PRIM_2); - } - - /*--- Limiter MPI ---*/ - - InitiateComms(geometry, config, PRIMITIVE_LIMITER); - CompleteComms(geometry, config, PRIMITIVE_LIMITER); - -#ifdef CODI_REVERSE_TYPE - if (TapeActive) AD::StartRecording(); -#endif + computeLimiters(kindLimiter, this, PRIMITIVE_LIMITER, PERIODIC_LIM_PRIM_1, PERIODIC_LIM_PRIM_2, + *geometry, *config, 0, nPrimVarGrad, primitives, gradient, primMin, primMax, limiter); } void CEulerSolver::SetPreconditioner(CConfig *config, unsigned long iPoint) { diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 12ff1babeee6..0e2a73da82ac 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -30,6 +30,10 @@ #include "../../Common/include/toolboxes/printing_toolbox.hpp" #include "../include/variables/CIncEulerVariable.hpp" #include "../include/variables/CIncNSVariable.hpp" +#include "../include/gradients/computeGradientsGreenGauss.hpp" +#include "../include/gradients/computeGradientsLeastSquares.hpp" +#include "../include/limiters/computeLimiters.hpp" + CIncEulerSolver::CIncEulerSolver(void) : CSolver() { /*--- Basic array initialization ---*/ @@ -3563,619 +3567,44 @@ void CIncEulerSolver::ComputeUnderRelaxationFactor(CSolver **solver_container, C } void CIncEulerSolver::SetPrimitive_Gradient_GG(CGeometry *geometry, CConfig *config, bool reconstruction) { - unsigned long iPoint, jPoint, iEdge, iVertex; - unsigned short iDim, iVar, iMarker; - su2double *PrimVar_Vertex, *PrimVar_i, *PrimVar_j, PrimVar_Average, - Partial_Gradient, Partial_Res, *Normal, Vol; - - /*--- Incompressible flow, primitive variables nDim+4, (P, vx, vy, vz, T, rho, beta) ---*/ - - PrimVar_Vertex = new su2double [nPrimVarGrad]; - PrimVar_i = new su2double [nPrimVarGrad]; - PrimVar_j = new su2double [nPrimVarGrad]; - - /*--- Set Gradient_Primitive to zero ---*/ - nodes->SetGradient_PrimitiveZero(); - - /*--- Loop interior edges ---*/ - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - PrimVar_i[iVar] = nodes->GetPrimitive(iPoint,iVar); - PrimVar_j[iVar] = nodes->GetPrimitive(jPoint,iVar); - } - - Normal = geometry->edge[iEdge]->GetNormal(); - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - PrimVar_Average = 0.5 * ( PrimVar_i[iVar] + PrimVar_j[iVar] ); - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Res = PrimVar_Average*Normal[iDim]; - if (geometry->node[iPoint]->GetDomain()) - nodes->AddGradient_Primitive(iPoint,iVar, iDim, Partial_Res); - if (geometry->node[jPoint]->GetDomain()) - nodes->SubtractGradient_Primitive(jPoint,iVar, iDim, Partial_Res); - } - } - } - - /*--- Loop boundary edges ---*/ - for (iMarker = 0; 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++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - if (geometry->node[iPoint]->GetDomain()) { - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - PrimVar_Vertex[iVar] = nodes->GetPrimitive(iPoint,iVar); - - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Res = PrimVar_Vertex[iVar]*Normal[iDim]; - nodes->SubtractGradient_Primitive(iPoint,iVar, iDim, Partial_Res); - } - } - } - } - } - - /*--- Correct the gradient values across any periodic boundaries. ---*/ - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_GG); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_GG); - } - - /*--- Update gradient value ---*/ - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Get the volume, which may include periodic components. ---*/ - - Vol = (geometry->node[iPoint]->GetVolume() + - geometry->node[iPoint]->GetPeriodicVolume()); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Gradient = nodes->GetGradient_Primitive(iPoint, iVar, iDim)/Vol; - if (reconstruction) - nodes->SetGradient_Reconstruction(iPoint, iVar, iDim, Partial_Gradient); - else - nodes->SetGradient_Primitive(iPoint, iVar, iDim, Partial_Gradient); - } - } - } - - /*--- Communicate the gradient values via MPI. ---*/ - - InitiateComms(geometry, config, PRIMITIVE_GRADIENT); - CompleteComms(geometry, config, PRIMITIVE_GRADIENT); - - delete [] PrimVar_Vertex; - delete [] PrimVar_i; - delete [] PrimVar_j; - + const auto& primitives = nodes->GetPrimitive(); + auto& gradient = reconstruction? nodes->GetGradient_Reconstruction() : nodes->GetGradient_Primitive(); + + computeGradientsGreenGauss(this, PRIMITIVE_GRADIENT, PERIODIC_PRIM_GG, *geometry, + *config, primitives, 0, nPrimVarGrad, gradient); } void CIncEulerSolver::SetPrimitive_Gradient_LS(CGeometry *geometry, CConfig *config, bool reconstruction) { - - unsigned short iVar, iDim, jDim, iNeigh; - unsigned long iPoint, jPoint; - su2double *PrimVar_i, *PrimVar_j, *Coord_i, *Coord_j; - su2double r11, r12, r13, r22, r23, r23_a, r23_b, r33, weight; - su2double z11, z12, z13, z22, z23, z33, detR2; - bool singular; - + /*--- Set a flag for unweighted or weighted least-squares. ---*/ - - bool weighted = true; - if (reconstruction) { - if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) - weighted = false; - } else if (config->GetKind_Gradient_Method() == LEAST_SQUARES) { - weighted = false; - } - - /*--- Clear Rmatrix, which could eventually be computed once - and stored for static meshes, as well as the prim gradient. ---*/ + bool weighted; - nodes->SetRmatrixZero(); - nodes->SetGradient_PrimitiveZero(); + if (reconstruction) + weighted = (config->GetKind_Gradient_Method_Recon() == WEIGHTED_LEAST_SQUARES); + else + weighted = (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES); - /*--- Loop over points of the grid ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Set the value of the singular ---*/ - singular = false; - - /*--- Get coordinates ---*/ - - Coord_i = geometry->node[iPoint]->GetCoord(); - - /*--- Get primitives from CVariable ---*/ - - PrimVar_i = nodes->GetPrimitive(iPoint); - - /*--- Inizialization of variables ---*/ - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) - for (iDim = 0; iDim < nDim; iDim++) - Cvector[iVar][iDim] = 0.0; - - /*--- Clear Rmatrix, which could eventually be computed once - and stored for static meshes, as well as the prim gradient. ---*/ - - AD::StartPreacc(); - AD::SetPreaccIn(PrimVar_i, nPrimVarGrad); - AD::SetPreaccIn(Coord_i, nDim); - - for (iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); iNeigh++) { - jPoint = geometry->node[iPoint]->GetPoint(iNeigh); - Coord_j = geometry->node[jPoint]->GetCoord(); - - PrimVar_j = nodes->GetPrimitive(jPoint); - - AD::SetPreaccIn(Coord_j, nDim); - AD::SetPreaccIn(PrimVar_j, nPrimVarGrad); - - if (weighted) { - weight = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - weight += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - } else { - weight = 1.0; - } - - /*--- Sumations for entries of upper triangular matrix R ---*/ - - if (weight != 0.0) { - - nodes->AddRmatrix(iPoint,0, 0, (Coord_j[0]-Coord_i[0])*(Coord_j[0]-Coord_i[0])/weight); - nodes->AddRmatrix(iPoint,0, 1, (Coord_j[0]-Coord_i[0])*(Coord_j[1]-Coord_i[1])/weight); - nodes->AddRmatrix(iPoint,1, 1, (Coord_j[1]-Coord_i[1])*(Coord_j[1]-Coord_i[1])/weight); - - if (nDim == 3) { - nodes->AddRmatrix(iPoint,0, 2, (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight); - nodes->AddRmatrix(iPoint,1, 2, (Coord_j[1]-Coord_i[1])*(Coord_j[2]-Coord_i[2])/weight); - nodes->AddRmatrix(iPoint,2, 1, (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight); - nodes->AddRmatrix(iPoint,2, 2, (Coord_j[2]-Coord_i[2])*(Coord_j[2]-Coord_i[2])/weight); - } - - /*--- Entries of c:= transpose(A)*b ---*/ - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - nodes->AddGradient_Primitive(iPoint,iVar,iDim, (Coord_j[iDim]-Coord_i[iDim])*(PrimVar_j[iVar]-PrimVar_i[iVar])/weight); - } - } - - } - } - AD::SetPreaccOut(nodes->GetRmatrix(iPoint), nDim, nDim); - AD::SetPreaccOut(nodes->GetGradient_Primitive(iPoint), nPrimVarGrad, nDim); - AD::EndPreacc(); - } - - /*--- Correct the gradient values across any periodic boundaries. ---*/ + const auto& primitives = nodes->GetPrimitive(); + auto& rmatrix = nodes->GetRmatrix(); + auto& gradient = reconstruction? nodes->GetGradient_Reconstruction() : nodes->GetGradient_Primitive(); + PERIODIC_QUANTITIES kindPeriodicComm = weighted? PERIODIC_PRIM_LS : PERIODIC_PRIM_ULS; - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - if (weighted) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_LS); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_LS); - } else { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_ULS); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_PRIM_ULS); - } - } - - /*--- Second loop over points of the grid to compute final gradient ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Set the value of the singular ---*/ - - singular = false; - - /*--- Entries of upper triangular matrix R ---*/ - - r11 = 0.0; r12 = 0.0; r13 = 0.0; r22 = 0.0; - r23 = 0.0; r23_a = 0.0; r23_b = 0.0; r33 = 0.0; - - r11 = nodes->GetRmatrix(iPoint,0,0); - r12 = nodes->GetRmatrix(iPoint,0,1); - r22 = nodes->GetRmatrix(iPoint,1,1); - - AD::StartPreacc(); - AD::SetPreaccIn(r11); - AD::SetPreaccIn(r12); - AD::SetPreaccIn(r22); - - if (r11 >= 0.0) r11 = sqrt(r11); else r11 = 0.0; - if (r11 != 0.0) r12 = r12/r11; else r12 = 0.0; - if (r22-r12*r12 >= 0.0) r22 = sqrt(r22-r12*r12); else r22 = 0.0; - - if (nDim == 3) { - r13 = nodes->GetRmatrix(iPoint,0,2); - r23_a = nodes->GetRmatrix(iPoint,1,2); - r23_b = nodes->GetRmatrix(iPoint,2,1); - r33 = nodes->GetRmatrix(iPoint,2,2); - - AD::SetPreaccIn(r13); - AD::SetPreaccIn(r23_a); - AD::SetPreaccIn(r23_b); - AD::SetPreaccIn(r33); - - if (r11 != 0.0) r13 = r13/r11; else r13 = 0.0; - if ((r22 != 0.0) && (r11*r22 != 0.0)) r23 = r23_a/r22 - r23_b*r12/(r11*r22); else r23 = 0.0; - if (r33-r23*r23-r13*r13 >= 0.0) r33 = sqrt(r33-r23*r23-r13*r13); else r33 = 0.0; - } - - /*--- Compute determinant ---*/ - - if (nDim == 2) detR2 = (r11*r22)*(r11*r22); - else detR2 = (r11*r22*r33)*(r11*r22*r33); - - /*--- Detect singular matrices ---*/ - - if (abs(detR2) <= EPS) { detR2 = 1.0; singular = true; } - - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - - if (singular) { - for (iDim = 0; iDim < nDim; iDim++) - for (jDim = 0; jDim < nDim; jDim++) - Smatrix[iDim][jDim] = 0.0; - } - else { - if (nDim == 2) { - Smatrix[0][0] = (r12*r12+r22*r22)/detR2; - Smatrix[0][1] = -r11*r12/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = r11*r11/detR2; - } - else { - z11 = r22*r33; z12 = -r12*r33; z13 = r12*r23-r13*r22; - z22 = r11*r33; z23 = -r11*r23; z33 = r11*r22; - Smatrix[0][0] = (z11*z11+z12*z12+z13*z13)/detR2; - Smatrix[0][1] = (z12*z22+z13*z23)/detR2; - Smatrix[0][2] = (z13*z33)/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = (z22*z22+z23*z23)/detR2; - Smatrix[1][2] = (z23*z33)/detR2; - Smatrix[2][0] = Smatrix[0][2]; - Smatrix[2][1] = Smatrix[1][2]; - Smatrix[2][2] = (z33*z33)/detR2; - } - } - - AD::SetPreaccOut(Smatrix, nDim, nDim); - AD::EndPreacc(); - - /*--- Computation of the gradient: S*c ---*/ - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - Cvector[iVar][iDim] = 0.0; - for (jDim = 0; jDim < nDim; jDim++) { - Cvector[iVar][iDim] += Smatrix[iDim][jDim]*nodes->GetGradient_Primitive(iPoint,iVar, jDim); - } - } - } - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - if (reconstruction) - nodes->SetGradient_Reconstruction(iPoint, iVar, iDim, Cvector[iVar][iDim]); - else - nodes->SetGradient_Primitive(iPoint, iVar, iDim, Cvector[iVar][iDim]); - } - } - - } - - /*--- Communicate the gradient values via MPI. ---*/ - - InitiateComms(geometry, config, PRIMITIVE_GRADIENT); - CompleteComms(geometry, config, PRIMITIVE_GRADIENT); - + computeGradientsLeastSquares(this, PRIMITIVE_GRADIENT, kindPeriodicComm, *geometry, *config, + weighted, primitives, 0, nPrimVarGrad, gradient, rmatrix); } void CIncEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) { - - unsigned long iEdge, iPoint, jPoint; - unsigned short iVar, iDim; - su2double **Gradient_i, **Gradient_j, *Coord_i, *Coord_j, - *Primitive, *Primitive_i, *Primitive_j, - *LocalMinPrimitive = NULL, *LocalMaxPrimitive = NULL, - *GlobalMinPrimitive = NULL, *GlobalMaxPrimitive = NULL, - dave, LimK, eps2, eps1, dm, dp, du, y, limiter; - -#ifdef CODI_REVERSE_TYPE - bool TapeActive = false; - - if (config->GetDiscrete_Adjoint() && config->GetFrozen_Limiter_Disc()) { - /*--- If limiters are frozen do not record the computation ---*/ - TapeActive = AD::globalTape.isActive(); - AD::StopRecording(); - } -#endif - - dave = config->GetRefElemLength(); - LimK = config->GetVenkat_LimiterCoeff(); - - if (config->GetKind_SlopeLimit_Flow() == NO_LIMITER) { - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - nodes->SetLimiter_Primitive(iPoint,iVar, 1.0); - } - } - - } - - else { - - /*--- Initialize solution max and solution min and the limiter in the entire domain --*/ - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - nodes->SetSolution_Max(iPoint,iVar, -EPS); - nodes->SetSolution_Min(iPoint,iVar, EPS); - nodes->SetLimiter_Primitive(iPoint,iVar, 2.0); - } - } - - /*--- Establish bounds for Spekreijse monotonicity by finding max & min values of neighbor variables --*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - - /*--- Get the primitive variables ---*/ - - Primitive_i = nodes->GetPrimitive(iPoint); - Primitive_j = nodes->GetPrimitive(jPoint); - - /*--- Compute the maximum, and minimum values for nodes i & j ---*/ - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - du = (Primitive_j[iVar] - Primitive_i[iVar]); - nodes->SetSolution_Min(iPoint, iVar, min(nodes->GetSolution_Min(iPoint, iVar), du)); - nodes->SetSolution_Max(iPoint, iVar, max(nodes->GetSolution_Max(iPoint, iVar), du)); - nodes->SetSolution_Min(jPoint, iVar, min(nodes->GetSolution_Min(jPoint, iVar), -du)); - nodes->SetSolution_Max(jPoint, iVar, max(nodes->GetSolution_Max(jPoint, iVar), -du)); - } - - } - - /*--- Correct the limiter values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_PRIM_1); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_PRIM_1); - } - - } - - - /*--- Barth-Jespersen limiter with Venkatakrishnan modification ---*/ - - if (config->GetKind_SlopeLimit_Flow() == BARTH_JESPERSEN) { - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Gradient_i = nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = nodes->GetGradient_Reconstruction(jPoint); - Coord_i = geometry->node[iPoint]->GetCoord(); - Coord_j = geometry->node[jPoint]->GetCoord(); - - AD::StartPreacc(); - AD::SetPreaccIn(Gradient_i, nPrimVarGrad, nDim); - AD::SetPreaccIn(Gradient_j, nPrimVarGrad, nDim); - AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - - AD::SetPreaccIn(nodes->GetSolution_Max(iPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Min(iPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Max(jPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Min(jPoint,iVar)); - - /*--- Calculate the interface left gradient, delta- (dm) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_j[iDim]-Coord_i[iDim])*Gradient_i[iVar][iDim]; - - if (dm == 0.0) { limiter = 2.0; } - else { - if ( dm > 0.0 ) dp = nodes->GetSolution_Max(iPoint,iVar); - else dp = nodes->GetSolution_Min(iPoint,iVar); - limiter = dp/dm; - } - - if (limiter < nodes->GetLimiter_Primitive(iPoint,iVar)) { - nodes->SetLimiter_Primitive(iPoint,iVar, limiter); - AD::SetPreaccOut(nodes->GetLimiter_Primitive(iPoint)[iVar]); - } - - /*--- Calculate the interface right gradient, delta+ (dp) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_i[iDim]-Coord_j[iDim])*Gradient_j[iVar][iDim]; - - if (dm == 0.0) { limiter = 2.0; } - else { - if ( dm > 0.0 ) dp = nodes->GetSolution_Max(jPoint,iVar); - else dp = nodes->GetSolution_Min(jPoint,iVar); - limiter = dp/dm; - } - - if (limiter < nodes->GetLimiter_Primitive(jPoint,iVar)) { - nodes->SetLimiter_Primitive(jPoint,iVar, limiter); - AD::SetPreaccOut(nodes->GetLimiter_Primitive(jPoint)[iVar]); - } - - } - - AD::EndPreacc(); - - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - y = nodes->GetLimiter_Primitive(iPoint,iVar); - limiter = (y*y + 2.0*y) / (y*y + y + 2.0); - nodes->SetLimiter_Primitive(iPoint,iVar, limiter); - } - } - - } - - /*--- Venkatakrishnan limiter ---*/ - - if ((config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN) || - (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG)) { - - if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { - - /*--- Allocate memory for the max and min primitive value --*/ - - LocalMinPrimitive = new su2double [nPrimVarGrad]; GlobalMinPrimitive = new su2double [nPrimVarGrad]; - LocalMaxPrimitive = new su2double [nPrimVarGrad]; GlobalMaxPrimitive = new su2double [nPrimVarGrad]; - - /*--- Compute the max value and min value of the solution ---*/ - - Primitive = nodes->GetPrimitive(0); - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - LocalMinPrimitive[iVar] = Primitive[iVar]; - LocalMaxPrimitive[iVar] = Primitive[iVar]; - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - - /*--- Get the primitive variables ---*/ - - Primitive = nodes->GetPrimitive(iPoint); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - LocalMinPrimitive[iVar] = min (LocalMinPrimitive[iVar], Primitive[iVar]); - LocalMaxPrimitive[iVar] = max (LocalMaxPrimitive[iVar], Primitive[iVar]); - } - - } - -#ifdef HAVE_MPI - SU2_MPI::Allreduce(LocalMinPrimitive, GlobalMinPrimitive, nPrimVarGrad, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - SU2_MPI::Allreduce(LocalMaxPrimitive, GlobalMaxPrimitive, nPrimVarGrad, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); -#else - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - GlobalMinPrimitive[iVar] = LocalMinPrimitive[iVar]; - GlobalMaxPrimitive[iVar] = LocalMaxPrimitive[iVar]; - } -#endif - } - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Gradient_i = nodes->GetGradient_Primitive(iPoint); - Gradient_j = nodes->GetGradient_Primitive(jPoint); - Coord_i = geometry->node[iPoint]->GetCoord(); - Coord_j = geometry->node[jPoint]->GetCoord(); - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - - AD::StartPreacc(); - AD::SetPreaccIn(Gradient_i[iVar], nDim); - AD::SetPreaccIn(Gradient_j[iVar], nDim); - AD::SetPreaccIn(Coord_i, nDim); - AD::SetPreaccIn(Coord_j, nDim); - AD::SetPreaccIn(nodes->GetSolution_Max(iPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Min(iPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Max(jPoint,iVar)); - AD::SetPreaccIn(nodes->GetSolution_Min(jPoint,iVar)); - - if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { - AD::SetPreaccIn(GlobalMaxPrimitive[iVar]); - AD::SetPreaccIn(GlobalMinPrimitive[iVar]); - eps1 = LimK * (GlobalMaxPrimitive[iVar] - GlobalMinPrimitive[iVar]); - eps2 = eps1*eps1; - } - else { - eps1 = LimK*dave; - eps2 = eps1*eps1*eps1; - } - - /*--- Calculate the interface left gradient, delta- (dm) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_j[iDim]-Coord_i[iDim])*Gradient_i[iVar][iDim]; - - /*--- Calculate the interface right gradient, delta+ (dp) ---*/ - - if ( dm > 0.0 ) dp = nodes->GetSolution_Max(iPoint,iVar); - else dp = nodes->GetSolution_Min(iPoint,iVar); - - limiter = ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < nodes->GetLimiter_Primitive(iPoint,iVar)){ - nodes->SetLimiter_Primitive(iPoint,iVar, limiter); - AD::SetPreaccOut(nodes->GetLimiter_Primitive(iPoint)[iVar]); - } - - /*-- Repeat for point j on the edge ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_i[iDim]-Coord_j[iDim])*Gradient_j[iVar][iDim]; - - if ( dm > 0.0 ) dp = nodes->GetSolution_Max(jPoint,iVar); - else dp = nodes->GetSolution_Min(jPoint,iVar); - - limiter = ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < nodes->GetLimiter_Primitive(jPoint,iVar)){ - nodes->SetLimiter_Primitive(jPoint,iVar, limiter); - AD::SetPreaccOut(nodes->GetLimiter_Primitive(jPoint)[iVar]); - } - - AD::EndPreacc(); - } - } - - if (LocalMinPrimitive != NULL) delete [] LocalMinPrimitive; - if (LocalMaxPrimitive != NULL) delete [] LocalMaxPrimitive; - if (GlobalMinPrimitive != NULL) delete [] GlobalMinPrimitive; - if (GlobalMaxPrimitive != NULL) delete [] GlobalMaxPrimitive; - } - - /*--- Correct the limiter values across any periodic boundaries. ---*/ + auto kindLimiter = static_cast(config->GetKind_SlopeLimit_Flow()); + const auto& primitives = nodes->GetPrimitive(); + const auto& gradient = nodes->GetGradient_Reconstruction(); + auto& primMin = nodes->GetSolution_Min(); + auto& primMax = nodes->GetSolution_Max(); + auto& limiter = nodes->GetLimiter_Primitive(); - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_PRIM_2); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_PRIM_2); - } - - /*--- Limiter MPI ---*/ - - InitiateComms(geometry, config, PRIMITIVE_LIMITER); - CompleteComms(geometry, config, PRIMITIVE_LIMITER); - -#ifdef CODI_REVERSE_TYPE - if (TapeActive) AD::StartRecording(); -#endif + computeLimiters(kindLimiter, this, PRIMITIVE_LIMITER, PERIODIC_LIM_PRIM_1, PERIODIC_LIM_PRIM_2, + *geometry, *config, 0, nPrimVarGrad, primitives, gradient, primMin, primMax, limiter); } void CIncEulerSolver::SetInletAtVertex(su2double *val_inlet, diff --git a/SU2_CFD/src/solver_structure.cpp b/SU2_CFD/src/solver_structure.cpp index d24189e10760..5ca64a2718cb 100644 --- a/SU2_CFD/src/solver_structure.cpp +++ b/SU2_CFD/src/solver_structure.cpp @@ -28,6 +28,9 @@ #include "../include/solver_structure.hpp" #include "../include/variables/CBaselineVariable.hpp" +#include "../include/gradients/computeGradientsGreenGauss.hpp" +#include "../include/gradients/computeGradientsLeastSquares.hpp" +#include "../include/limiters/computeLimiters.hpp" #include "../../Common/include/toolboxes/MMS/CIncTGVSolution.hpp" #include "../../Common/include/toolboxes/MMS/CInviscidVortexSolution.hpp" #include "../../Common/include/toolboxes/MMS/CMMSIncEulerSolution.hpp" @@ -267,9 +270,13 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, CConfig *config, unsigned short val_periodic_index, unsigned short commType) { - + + /*--- Check for dummy communication. ---*/ + + if (commType == PERIODIC_NONE) return; + /*--- Local variables ---*/ - + bool boundary_i, boundary_j; bool weighted = true; @@ -1284,12 +1291,24 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, ensures that the proper min and max of the solution are found among all nodes adjacent to periodic faces. ---*/ + /*--- We send the min and max over "our" neighbours. ---*/ + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { Sol_Min[iVar] = base_nodes->GetSolution_Min(iPoint, iVar); Sol_Max[iVar] = base_nodes->GetSolution_Max(iPoint, iVar); - - bufDSend[buf_offset+iVar] = base_nodes->GetSolution_Min(iPoint, iVar); - bufDSend[buf_offset+nPrimVarGrad+iVar] = base_nodes->GetSolution_Max(iPoint, iVar); + } + + for (iNeighbor = 0; iNeighbor < geometry->node[iPoint]->GetnPoint(); iNeighbor++) { + jPoint = geometry->node[iPoint]->GetPoint(iNeighbor); + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + Sol_Min[iVar] = min(Sol_Min[iVar], base_nodes->GetPrimitive(jPoint, iVar)); + Sol_Max[iVar] = max(Sol_Max[iVar], base_nodes->GetPrimitive(jPoint, iVar)); + } + } + + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + bufDSend[buf_offset+iVar] = Sol_Min[iVar]; + bufDSend[buf_offset+nPrimVarGrad+iVar] = Sol_Max[iVar]; } /*--- Rotate the momentum components of the min/max. ---*/ @@ -1370,12 +1389,24 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, ensures that the proper min and max of the solution are found among all nodes adjacent to periodic faces. ---*/ + /*--- We send the min and max over "our" neighbours. ---*/ + for (iVar = 0; iVar < nVar; iVar++) { Sol_Min[iVar] = base_nodes->GetSolution_Min(iPoint, iVar); Sol_Max[iVar] = base_nodes->GetSolution_Max(iPoint, iVar); - - bufDSend[buf_offset+iVar] = base_nodes->GetSolution_Min(iPoint, iVar); - bufDSend[buf_offset+nVar+iVar] = base_nodes->GetSolution_Max(iPoint, iVar); + } + + for (iNeighbor = 0; iNeighbor < geometry->node[iPoint]->GetnPoint(); iNeighbor++) { + jPoint = geometry->node[iPoint]->GetPoint(iNeighbor); + for (iVar = 0; iVar < nVar; iVar++) { + Sol_Min[iVar] = min(Sol_Min[iVar], base_nodes->GetSolution(jPoint, iVar)); + Sol_Max[iVar] = max(Sol_Max[iVar], base_nodes->GetSolution(jPoint, iVar)); + } + } + + for (iVar = 0; iVar < nVar; iVar++) { + bufDSend[buf_offset+iVar] = Sol_Min[iVar]; + bufDSend[buf_offset+nVar+iVar] = Sol_Max[iVar]; } /*--- Rotate the momentum components of the min/max. ---*/ @@ -1487,9 +1518,13 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, CConfig *config, unsigned short val_periodic_index, unsigned short commType) { - + + /*--- Check for dummy communication. ---*/ + + if (commType == PERIODIC_NONE) return; + /*--- Local variables ---*/ - + unsigned short nPeriodic = config->GetnMarker_Periodic(); unsigned short iDim, jDim, iVar, jVar, iPeriodic, nNeighbor; @@ -1769,13 +1804,23 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, case PERIODIC_LIM_PRIM_1: - /*--- Check the min and max values found on the matching - perioic faces for the solution, and store the proper min - and max for this point. ---*/ + /*--- Update solution min/max with min/max between "us" and + the periodic match plus its neighbors, computation will need to + be concluded on "our" side to account for "our" neighbors. ---*/ for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - base_nodes->SetSolution_Min(iPoint, iVar, min(base_nodes->GetSolution_Min(iPoint, iVar), bufDRecv[buf_offset+iVar])); - base_nodes->SetSolution_Max(iPoint, iVar, max(base_nodes->GetSolution_Max(iPoint, iVar), bufDRecv[buf_offset+nPrimVarGrad+iVar])); + + /*--- Solution minimum. ---*/ + + Solution_Min = min(base_nodes->GetSolution_Min(iPoint, iVar), + bufDRecv[buf_offset+iVar]); + base_nodes->SetSolution_Min(iPoint, iVar, Solution_Min); + + /*--- Solution maximum. ---*/ + + Solution_Max = max(base_nodes->GetSolution_Max(iPoint, iVar), + bufDRecv[buf_offset+nPrimVarGrad+iVar]); + base_nodes->SetSolution_Max(iPoint, iVar, Solution_Max); } break; @@ -1786,16 +1831,18 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, faces for the limiter, and store the proper min value. ---*/ for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - base_nodes->SetLimiter_Primitive(iPoint, iVar, min(base_nodes->GetLimiter_Primitive(iPoint, iVar), bufDRecv[buf_offset+iVar])); + Limiter_Min = min(base_nodes->GetLimiter_Primitive(iPoint, iVar), + bufDRecv[buf_offset+iVar]); + base_nodes->SetLimiter_Primitive(iPoint, iVar, Limiter_Min); } break; case PERIODIC_LIM_SOL_1: - /*--- Check the min and max values found on the matching - perioic faces for the solution, and store the proper min - and max for this point. ---*/ + /*--- Update solution min/max with min/max between "us" and + the periodic match plus its neighbors, computation will need to + be concluded on "our" side to account for "our" neighbors. ---*/ for (iVar = 0; iVar < nVar; iVar++) { @@ -2754,472 +2801,51 @@ void CSolver::SetRotatingFrame_GCL(CGeometry *geometry, CConfig *config) { } void CSolver::SetAuxVar_Gradient_GG(CGeometry *geometry, CConfig *config) { - - unsigned long Point = 0, iPoint = 0, jPoint = 0, iEdge, iVertex; - unsigned short nDim = geometry->GetnDim(), iDim, iMarker; - - su2double AuxVar_Vertex, AuxVar_i, AuxVar_j, AuxVar_Average; - su2double *Gradient, DualArea, Partial_Res, Grad_Val, *Normal; - - /*--- Set Gradient to Zero ---*/ - base_nodes->SetAuxVarGradientZero(); - - /*--- Loop interior edges ---*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - - AuxVar_i = base_nodes->GetAuxVar(iPoint); - AuxVar_j = base_nodes->GetAuxVar(jPoint); - - Normal = geometry->edge[iEdge]->GetNormal(); - AuxVar_Average = 0.5 * ( AuxVar_i + AuxVar_j); - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Res = AuxVar_Average*Normal[iDim]; - base_nodes->AddAuxVarGradient(iPoint, iDim, Partial_Res); - base_nodes->SubtractAuxVarGradient(jPoint,iDim, Partial_Res); - } - } - - /*--- Loop boundary edges ---*/ - - for (iMarker = 0; 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++) { - Point = geometry->vertex[iMarker][iVertex]->GetNode(); - AuxVar_Vertex = base_nodes->GetAuxVar(Point); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Res = AuxVar_Vertex*Normal[iDim]; - base_nodes->SubtractAuxVarGradient(Point,iDim, Partial_Res); - } - } - } - - for (iPoint=0; iPointGetnPoint(); iPoint++) - for (iDim = 0; iDim < nDim; iDim++) { - Gradient = base_nodes->GetAuxVarGradient(iPoint); - DualArea = geometry->node[iPoint]->GetVolume(); - Grad_Val = Gradient[iDim]/(DualArea+EPS); - base_nodes->SetAuxVarGradient(iPoint, iDim, Grad_Val); - } - - /*--- Gradient MPI ---*/ - - InitiateComms(geometry, config, AUXVAR_GRADIENT); - CompleteComms(geometry, config, AUXVAR_GRADIENT); + const auto solution = base_nodes->GetAuxVar(); + auto gradient = base_nodes->GetAuxVarGradient(); + computeGradientsGreenGauss(this, AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, + *config, solution, 0, 1, gradient); } void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, CConfig *config) { - - unsigned short iDim, jDim, iNeigh; - unsigned short nDim = geometry->GetnDim(); - unsigned long iPoint, jPoint; - su2double *Coord_i, *Coord_j, AuxVar_i, AuxVar_j, weight, r11, r12, r13, r22, r23, r23_a, - r23_b, r33, z11, z12, z13, z22, z23, z33, detR2, product; - bool singular = false; - - su2double *Cvector = new su2double [nDim]; - - /*--- Loop over points of the grid ---*/ - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - - Coord_i = geometry->node[iPoint]->GetCoord(); - AuxVar_i = base_nodes->GetAuxVar(iPoint); - - /*--- Inizialization of variables ---*/ - for (iDim = 0; iDim < nDim; iDim++) - Cvector[iDim] = 0.0; - - r11 = 0.0; r12 = 0.0; r13 = 0.0; r22 = 0.0; - r23 = 0.0; r23_a = 0.0; r23_b = 0.0; r33 = 0.0; - - for (iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); iNeigh++) { - jPoint = geometry->node[iPoint]->GetPoint(iNeigh); - Coord_j = geometry->node[jPoint]->GetCoord(); - AuxVar_j = base_nodes->GetAuxVar(jPoint); - - weight = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - weight += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - - /*--- Sumations for entries of upper triangular matrix R ---*/ - - if (fabs(weight) > EPS) { - r11 += (Coord_j[0]-Coord_i[0])*(Coord_j[0]-Coord_i[0])/weight; - r12 += (Coord_j[0]-Coord_i[0])*(Coord_j[1]-Coord_i[1])/weight; - r22 += (Coord_j[1]-Coord_i[1])*(Coord_j[1]-Coord_i[1])/weight; - if (nDim == 3) { - r13 += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight; - r23_a += (Coord_j[1]-Coord_i[1])*(Coord_j[2]-Coord_i[2])/weight; - r23_b += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight; - r33 += (Coord_j[2]-Coord_i[2])*(Coord_j[2]-Coord_i[2])/weight; - } - - /*--- Entries of c:= transpose(A)*b ---*/ - - for (iDim = 0; iDim < nDim; iDim++) - Cvector[iDim] += (Coord_j[iDim]-Coord_i[iDim])*(AuxVar_j-AuxVar_i)/(weight); - } - - } - - /*--- Entries of upper triangular matrix R ---*/ - - if (fabs(r11) < EPS) r11 = EPS; - r11 = sqrt(r11); - r12 = r12/r11; - r22 = sqrt(r22-r12*r12); - if (fabs(r22) < EPS) r22 = EPS; - if (nDim == 3) { - r13 = r13/r11; - r23 = r23_a/(r22) - r23_b*r12/(r11*r22); - r33 = sqrt(r33-r23*r23-r13*r13); - } - - /*--- Compute determinant ---*/ - - if (nDim == 2) detR2 = (r11*r22)*(r11*r22); - else detR2 = (r11*r22*r33)*(r11*r22*r33); - - /*--- Detect singular matrices ---*/ - - if (fabs(detR2) < EPS) singular = true; - - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - - if (singular) { - for (iDim = 0; iDim < nDim; iDim++) - for (jDim = 0; jDim < nDim; jDim++) - Smatrix[iDim][jDim] = 0.0; - } - else { - if (nDim == 2) { - Smatrix[0][0] = (r12*r12+r22*r22)/detR2; - Smatrix[0][1] = -r11*r12/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = r11*r11/detR2; - } - else { - z11 = r22*r33; z12 = -r12*r33; z13 = r12*r23-r13*r22; - z22 = r11*r33; z23 = -r11*r23; z33 = r11*r22; - Smatrix[0][0] = (z11*z11+z12*z12+z13*z13)/detR2; - Smatrix[0][1] = (z12*z22+z13*z23)/detR2; - Smatrix[0][2] = (z13*z33)/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = (z22*z22+z23*z23)/detR2; - Smatrix[1][2] = (z23*z33)/detR2; - Smatrix[2][0] = Smatrix[0][2]; - Smatrix[2][1] = Smatrix[1][2]; - Smatrix[2][2] = (z33*z33)/detR2; - } - } - - /*--- Computation of the gradient: S*c ---*/ - - for (iDim = 0; iDim < nDim; iDim++) { - product = 0.0; - for (jDim = 0; jDim < nDim; jDim++) - product += Smatrix[iDim][jDim]*Cvector[jDim]; - if (geometry->node[iPoint]->GetDomain()) - base_nodes->SetAuxVarGradient(iPoint, iDim, product); - } - } - - delete [] Cvector; - - /*--- Gradient MPI ---*/ - - InitiateComms(geometry, config, AUXVAR_GRADIENT); - CompleteComms(geometry, config, AUXVAR_GRADIENT); - + + bool weighted = true; + const auto solution = base_nodes->GetAuxVar(); + auto gradient = base_nodes->GetAuxVarGradient(); + auto& rmatrix = base_nodes->GetRmatrix(); + + computeGradientsLeastSquares(this, AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config, + weighted, solution, 0, 1, gradient, rmatrix); } void CSolver::SetSolution_Gradient_GG(CGeometry *geometry, CConfig *config, bool reconstruction) { - unsigned long Point = 0, iPoint = 0, jPoint = 0, iEdge, iVertex; - unsigned short iVar, iDim, iMarker; - su2double *Solution_Vertex, *Solution_i, *Solution_j, Solution_Average, **Gradient, - Partial_Res, Grad_Val, *Normal, Vol; - - /*--- Set Gradient to Zero ---*/ - base_nodes->SetGradientZero(); - - /*--- Loop interior edges ---*/ - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - - Solution_i = base_nodes->GetSolution(iPoint); - Solution_j = base_nodes->GetSolution(jPoint); - Normal = geometry->edge[iEdge]->GetNormal(); - for (iVar = 0; iVar< nVar; iVar++) { - Solution_Average = 0.5 * (Solution_i[iVar] + Solution_j[iVar]); - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Res = Solution_Average*Normal[iDim]; - if (geometry->node[iPoint]->GetDomain()) - base_nodes->AddGradient(iPoint, iVar, iDim, Partial_Res); - if (geometry->node[jPoint]->GetDomain()) - base_nodes->SubtractGradient(jPoint,iVar, iDim, Partial_Res); - } - } - } - - /*--- Loop boundary edges ---*/ - for (iMarker = 0; 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++) { - Point = geometry->vertex[iMarker][iVertex]->GetNode(); - Solution_Vertex = base_nodes->GetSolution(Point); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - for (iVar = 0; iVar < nVar; iVar++) - for (iDim = 0; iDim < nDim; iDim++) { - Partial_Res = Solution_Vertex[iVar]*Normal[iDim]; - if (geometry->node[Point]->GetDomain()) - base_nodes->SubtractGradient(Point,iVar, iDim, Partial_Res); - } - } - } - } - - /*--- Correct the gradient values for any periodic boundaries. ---*/ - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_SOL_GG); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_SOL_GG); - } - - /*--- Compute gradient ---*/ - for (iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++) { - - /*--- Get the volume, which may include periodic components. ---*/ - - Vol = (geometry->node[iPoint]->GetVolume() + - geometry->node[iPoint]->GetPeriodicVolume()); - - for (iVar = 0; iVar < nVar; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - Gradient = base_nodes->GetGradient(iPoint); - Grad_Val = Gradient[iVar][iDim] / (Vol+EPS); - if (reconstruction) - base_nodes->SetGradient_Reconstruction(iPoint, iVar, iDim, Grad_Val); - else - base_nodes->SetGradient(iPoint, iVar, iDim, Grad_Val); - } - } - - } - - /*--- Gradient MPI ---*/ - - InitiateComms(geometry, config, SOLUTION_GRADIENT); - CompleteComms(geometry, config, SOLUTION_GRADIENT); - + const auto& solution = base_nodes->GetSolution(); + auto& gradient = reconstruction? base_nodes->GetGradient_Reconstruction() : base_nodes->GetGradient(); + + computeGradientsGreenGauss(this, SOLUTION_GRADIENT, PERIODIC_SOL_GG, *geometry, + *config, solution, 0, nVar, gradient); } void CSolver::SetSolution_Gradient_LS(CGeometry *geometry, CConfig *config, bool reconstruction) { - - unsigned short iDim, jDim, iVar, iNeigh; - unsigned long iPoint, jPoint; - su2double *Coord_i, *Coord_j, *Solution_i, *Solution_j; - su2double r11, r12, r13, r22, r23, r23_a, r23_b, r33, weight; - su2double detR2, z11, z12, z13, z22, z23, z33; - bool singular = false; /*--- Set a flag for unweighted or weighted least-squares. ---*/ - - bool weighted = true; - if (reconstruction) { - if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) - weighted = false; - } else if (config->GetKind_Gradient_Method() == LEAST_SQUARES) { - weighted = false; - } - - /*--- Clear Rmatrix, which could eventually be computed once - and stored for static meshes, as well as the gradient. ---*/ + bool weighted; - base_nodes->SetRmatrixZero(); - base_nodes->SetGradientZero(); + if (reconstruction) + weighted = (config->GetKind_Gradient_Method_Recon() == WEIGHTED_LEAST_SQUARES); + else + weighted = (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES); - /*--- Loop over points of the grid ---*/ - - for (iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++) { - - /*--- Set the value of the singular ---*/ - singular = false; - - /*--- Get coordinates ---*/ - - Coord_i = geometry->node[iPoint]->GetCoord(); - - /*--- Get consevative solution ---*/ - - Solution_i = base_nodes->GetSolution(iPoint); - - /*--- Inizialization of variables ---*/ - - for (iVar = 0; iVar < nVar; iVar++) - for (iDim = 0; iDim < nDim; iDim++) - Cvector[iVar][iDim] = 0.0; - - for (iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); iNeigh++) { - jPoint = geometry->node[iPoint]->GetPoint(iNeigh); - Coord_j = geometry->node[jPoint]->GetCoord(); - - Solution_j = base_nodes->GetSolution(jPoint); + const auto& solution = base_nodes->GetSolution(); + auto& rmatrix = base_nodes->GetRmatrix(); + auto& gradient = reconstruction? base_nodes->GetGradient_Reconstruction() : base_nodes->GetGradient(); + PERIODIC_QUANTITIES kindPeriodicComm = weighted? PERIODIC_SOL_LS : PERIODIC_SOL_ULS; - if (weighted) { - weight = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - weight += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - } else { - weight = 1.0; - } - - /*--- Sumations for entries of upper triangular matrix R ---*/ - - if (weight != 0.0) { - - base_nodes->AddRmatrix(iPoint,0, 0, (Coord_j[0]-Coord_i[0])*(Coord_j[0]-Coord_i[0])/weight); - base_nodes->AddRmatrix(iPoint,0, 1, (Coord_j[0]-Coord_i[0])*(Coord_j[1]-Coord_i[1])/weight); - base_nodes->AddRmatrix(iPoint,1, 1, (Coord_j[1]-Coord_i[1])*(Coord_j[1]-Coord_i[1])/weight); - - if (nDim == 3) { - base_nodes->AddRmatrix(iPoint,0, 2, (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight); - base_nodes->AddRmatrix(iPoint,1, 2, (Coord_j[1]-Coord_i[1])*(Coord_j[2]-Coord_i[2])/weight); - base_nodes->AddRmatrix(iPoint,2, 1, (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight); - base_nodes->AddRmatrix(iPoint,2, 2, (Coord_j[2]-Coord_i[2])*(Coord_j[2]-Coord_i[2])/weight); - } - - /*--- Entries of c:= transpose(A)*b ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - base_nodes->AddGradient(iPoint, iVar,iDim, (Coord_j[iDim]-Coord_i[iDim])*(Solution_j[iVar]-Solution_i[iVar])/weight); - } - } - - } - } - } - - /*--- Correct the gradient values for any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - if (weighted) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_SOL_LS); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_SOL_LS); - } else { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_SOL_ULS); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_SOL_ULS); - } - } - - /*--- Second loop over points of the grid to compute final gradient ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Set the value of the singular ---*/ - - singular = false; - - /*--- Entries of upper triangular matrix R ---*/ - - r11 = 0.0; r12 = 0.0; r13 = 0.0; r22 = 0.0; - r23 = 0.0; r23_a = 0.0; r23_b = 0.0; r33 = 0.0; - - r11 = base_nodes->GetRmatrix(iPoint,0,0); - r12 = base_nodes->GetRmatrix(iPoint,0,1); - r22 = base_nodes->GetRmatrix(iPoint,1,1); - - /*--- Entries of upper triangular matrix R ---*/ - - if (r11 >= 0.0) r11 = sqrt(r11); else r11 = 0.0; - if (r11 != 0.0) r12 = r12/r11; else r12 = 0.0; - if (r22-r12*r12 >= 0.0) r22 = sqrt(r22-r12*r12); else r22 = 0.0; - - if (nDim == 3) { - r13 = base_nodes->GetRmatrix(iPoint,0,2); - r23_a = base_nodes->GetRmatrix(iPoint,1,2); - r23_b = base_nodes->GetRmatrix(iPoint,2,1); - r33 = base_nodes->GetRmatrix(iPoint,2,2); - - if (r11 != 0.0) r13 = r13/r11; else r13 = 0.0; - if ((r22 != 0.0) && (r11*r22 != 0.0)) r23 = r23_a/r22 - r23_b*r12/(r11*r22); else r23 = 0.0; - if (r33-r23*r23-r13*r13 >= 0.0) r33 = sqrt(r33-r23*r23-r13*r13); else r33 = 0.0; - } - - /*--- Compute determinant ---*/ - - if (nDim == 2) detR2 = (r11*r22)*(r11*r22); - else detR2 = (r11*r22*r33)*(r11*r22*r33); - - /*--- Detect singular matrices ---*/ - - if (abs(detR2) <= EPS) { detR2 = 1.0; singular = true; } - - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - - if (singular) { - for (iDim = 0; iDim < nDim; iDim++) - for (jDim = 0; jDim < nDim; jDim++) - Smatrix[iDim][jDim] = 0.0; - } - else { - if (nDim == 2) { - Smatrix[0][0] = (r12*r12+r22*r22)/detR2; - Smatrix[0][1] = -r11*r12/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = r11*r11/detR2; - } - else { - z11 = r22*r33; z12 = -r12*r33; z13 = r12*r23-r13*r22; - z22 = r11*r33; z23 = -r11*r23; z33 = r11*r22; - Smatrix[0][0] = (z11*z11+z12*z12+z13*z13)/detR2; - Smatrix[0][1] = (z12*z22+z13*z23)/detR2; - Smatrix[0][2] = (z13*z33)/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = (z22*z22+z23*z23)/detR2; - Smatrix[1][2] = (z23*z33)/detR2; - Smatrix[2][0] = Smatrix[0][2]; - Smatrix[2][1] = Smatrix[1][2]; - Smatrix[2][2] = (z33*z33)/detR2; - } - } - - /*--- Computation of the gradient: S*c ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - Cvector[iVar][iDim] = 0.0; - for (jDim = 0; jDim < nDim; jDim++) { - Cvector[iVar][iDim] += Smatrix[iDim][jDim]*base_nodes->GetGradient(iPoint, iVar, jDim); - } - } - } - - for (iVar = 0; iVar < nVar; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { - if (reconstruction) - base_nodes->SetGradient_Reconstruction(iPoint, iVar, iDim, Cvector[iVar][iDim]); - else - base_nodes->SetGradient(iPoint, iVar, iDim, Cvector[iVar][iDim]); - } - } - - } - - /*--- Gradient MPI ---*/ - - InitiateComms(geometry, config, SOLUTION_GRADIENT); - CompleteComms(geometry, config, SOLUTION_GRADIENT); - + computeGradientsLeastSquares(this, SOLUTION_GRADIENT, kindPeriodicComm, *geometry, *config, + weighted, solution, 0, nVar, gradient, rmatrix); } void CSolver::Add_External_To_Solution() { @@ -3483,447 +3109,16 @@ void CSolver::SetAuxVar_Surface_Gradient(CGeometry *geometry, CConfig *config) { } void CSolver::SetSolution_Limiter(CGeometry *geometry, CConfig *config) { - - unsigned long iEdge, iPoint, jPoint; - unsigned short iVar, iDim; - su2double **Gradient_i, **Gradient_j, *Coord_i, *Coord_j, - *Solution, *Solution_i, *Solution_j, - *LocalMinSolution = NULL, *LocalMaxSolution = NULL, - *GlobalMinSolution = NULL, *GlobalMaxSolution = NULL, - dave, LimK, eps1, eps2, dm, dp, du, ds, y, limiter, SharpEdge_Distance; - -#ifdef CODI_REVERSE_TYPE - bool TapeActive = false; - - if (config->GetDiscrete_Adjoint() && config->GetFrozen_Limiter_Disc()) { - /*--- If limiters are frozen do not record the computation ---*/ - TapeActive = AD::globalTape.isActive(); - AD::StopRecording(); - } -#endif - - dave = config->GetRefElemLength(); - LimK = config->GetVenkat_LimiterCoeff(); - - if (config->GetKind_SlopeLimit() == NO_LIMITER) { - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - base_nodes->SetLimiter(iPoint, iVar, 1.0); - } - } - - } - - else { - - /*--- Initialize solution max and solution min and the limiter in the entire domain --*/ - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - base_nodes->SetSolution_Max(iPoint, iVar, -EPS); - base_nodes->SetSolution_Min(iPoint, iVar, EPS); - base_nodes->SetLimiter(iPoint, iVar, 2.0); - } - } - - /*--- Establish bounds for Spekreijse monotonicity by finding max & min values of neighbor variables --*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - - /*--- Get the conserved variables ---*/ - - Solution_i = base_nodes->GetSolution(iPoint); - Solution_j = base_nodes->GetSolution(jPoint); - - /*--- Compute the maximum, and minimum values for nodes i & j ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - du = (Solution_j[iVar] - Solution_i[iVar]); - base_nodes->SetSolution_Min(iPoint, iVar, min(base_nodes->GetSolution_Min(iPoint, iVar), du)); - base_nodes->SetSolution_Max(iPoint, iVar, max(base_nodes->GetSolution_Max(iPoint, iVar), du)); - base_nodes->SetSolution_Min(jPoint, iVar, min(base_nodes->GetSolution_Min(jPoint, iVar), -du)); - base_nodes->SetSolution_Max(jPoint, iVar, max(base_nodes->GetSolution_Max(jPoint, iVar), -du)); - } - - } - - /*--- Correct the limiter values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_SOL_1); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_SOL_1); - } - - } - - /*--- Barth-Jespersen limiter with Venkatakrishnan modification ---*/ - - if (config->GetKind_SlopeLimit_Flow() == BARTH_JESPERSEN) { - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Gradient_i = base_nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = base_nodes->GetGradient_Reconstruction(jPoint); - Coord_i = geometry->node[iPoint]->GetCoord(); - Coord_j = geometry->node[jPoint]->GetCoord(); - - AD::StartPreacc(); - AD::SetPreaccIn(Gradient_i, nVar, nDim); - AD::SetPreaccIn(Gradient_j, nVar, nDim); - AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); - - for (iVar = 0; iVar < nVar; iVar++) { - - AD::SetPreaccIn(base_nodes->GetSolution_Max(iPoint, iVar)); - AD::SetPreaccIn(base_nodes->GetSolution_Min(iPoint, iVar)); - AD::SetPreaccIn(base_nodes->GetSolution_Max(jPoint,iVar)); - AD::SetPreaccIn(base_nodes->GetSolution_Min(jPoint,iVar)); - /*--- Calculate the interface left gradient, delta- (dm) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_j[iDim]-Coord_i[iDim])*Gradient_i[iVar][iDim]; - - if (dm == 0.0) { limiter = 2.0; } - else { - if ( dm > 0.0 ) dp = base_nodes->GetSolution_Max(iPoint, iVar); - else dp = base_nodes->GetSolution_Min(iPoint, iVar); - limiter = dp/dm; - } - - if (limiter < base_nodes->GetLimiter(iPoint, iVar)) { - base_nodes->SetLimiter(iPoint, iVar, limiter); - AD::SetPreaccOut(base_nodes->GetLimiter(iPoint)[iVar]); - } - - /*--- Calculate the interface right gradient, delta+ (dp) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_i[iDim]-Coord_j[iDim])*Gradient_j[iVar][iDim]; - - if (dm == 0.0) { limiter = 2.0; } - else { - if ( dm > 0.0 ) dp = base_nodes->GetSolution_Max(jPoint,iVar); - else dp = base_nodes->GetSolution_Min(jPoint,iVar); - limiter = dp/dm; - } - - if (limiter < base_nodes->GetLimiter(jPoint,iVar)) { - base_nodes->SetLimiter(jPoint,iVar, limiter); - AD::SetPreaccOut(base_nodes->GetLimiter(jPoint)[iVar]); - } - - } - - AD::EndPreacc(); - - } - - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - y = base_nodes->GetLimiter(iPoint, iVar); - limiter = (y*y + 2.0*y) / (y*y + y + 2.0); - base_nodes->SetLimiter(iPoint, iVar, limiter); - } - } - - } - - /*--- Venkatakrishnan limiter ---*/ - - if ((config->GetKind_SlopeLimit() == VENKATAKRISHNAN) || (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG)) { - - if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { + auto kindLimiter = static_cast(config->GetKind_SlopeLimit()); + const auto& solution = base_nodes->GetSolution(); + const auto& gradient = base_nodes->GetGradient_Reconstruction(); + auto& solMin = base_nodes->GetSolution_Min(); + auto& solMax = base_nodes->GetSolution_Max(); + auto& limiter = base_nodes->GetLimiter(); - /*--- Allocate memory for the max and min solution value --*/ - - LocalMinSolution = new su2double [nVar]; GlobalMinSolution = new su2double [nVar]; - LocalMaxSolution = new su2double [nVar]; GlobalMaxSolution = new su2double [nVar]; - - /*--- Compute the max value and min value of the solution ---*/ - - Solution = base_nodes->GetSolution(iPoint); - for (iVar = 0; iVar < nVar; iVar++) { - LocalMinSolution[iVar] = Solution[iVar]; - LocalMaxSolution[iVar] = Solution[iVar]; - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - - /*--- Get the solution variables ---*/ - - Solution = base_nodes->GetSolution(iPoint); - - for (iVar = 0; iVar < nVar; iVar++) { - LocalMinSolution[iVar] = min (LocalMinSolution[iVar], Solution[iVar]); - LocalMaxSolution[iVar] = max (LocalMaxSolution[iVar], Solution[iVar]); - } - - } - -#ifdef HAVE_MPI - SU2_MPI::Allreduce(LocalMinSolution, GlobalMinSolution, nVar, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - SU2_MPI::Allreduce(LocalMaxSolution, GlobalMaxSolution, nVar, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); -#else - for (iVar = 0; iVar < nVar; iVar++) { - GlobalMinSolution[iVar] = LocalMinSolution[iVar]; - GlobalMaxSolution[iVar] = LocalMaxSolution[iVar]; - } -#endif - } - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Gradient_i = base_nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = base_nodes->GetGradient_Reconstruction(jPoint); - Coord_i = geometry->node[iPoint]->GetCoord(); - Coord_j = geometry->node[jPoint]->GetCoord(); - - AD::StartPreacc(); - AD::SetPreaccIn(Gradient_i, nVar, nDim); - AD::SetPreaccIn(Gradient_j, nVar, nDim); - AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); - - for (iVar = 0; iVar < nVar; iVar++) { - - AD::StartPreacc(); - AD::SetPreaccIn(Gradient_i[iVar], nDim); - AD::SetPreaccIn(Gradient_j[iVar], nDim); - AD::SetPreaccIn(Coord_i, nDim); - AD::SetPreaccIn(Coord_j, nDim); - AD::SetPreaccIn(base_nodes->GetSolution_Max(iPoint, iVar)); - AD::SetPreaccIn(base_nodes->GetSolution_Min(iPoint, iVar)); - AD::SetPreaccIn(base_nodes->GetSolution_Max(jPoint,iVar)); - AD::SetPreaccIn(base_nodes->GetSolution_Min(jPoint,iVar)); - - if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { - AD::SetPreaccIn(GlobalMaxSolution[iVar]); - AD::SetPreaccIn(GlobalMinSolution[iVar]); - eps1 = LimK * (GlobalMaxSolution[iVar] - GlobalMinSolution[iVar]); - eps2 = eps1*eps1; - } - else { - eps1 = LimK*dave; - eps2 = eps1*eps1*eps1; - } - - /*--- Calculate the interface left gradient, delta- (dm) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_j[iDim]-Coord_i[iDim])*Gradient_i[iVar][iDim]; - - /*--- Calculate the interface right gradient, delta+ (dp) ---*/ - - if ( dm > 0.0 ) dp = base_nodes->GetSolution_Max(iPoint, iVar); - else dp = base_nodes->GetSolution_Min(iPoint, iVar); - - limiter = ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < base_nodes->GetLimiter(iPoint, iVar)) { - base_nodes->SetLimiter(iPoint, iVar, limiter); - AD::SetPreaccOut(base_nodes->GetLimiter(iPoint)[iVar]); - } - - /*-- Repeat for point j on the edge ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_i[iDim]-Coord_j[iDim])*Gradient_j[iVar][iDim]; - - if ( dm > 0.0 ) dp = base_nodes->GetSolution_Max(jPoint,iVar); - else dp = base_nodes->GetSolution_Min(jPoint,iVar); - - limiter = ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < base_nodes->GetLimiter(jPoint,iVar)) { - base_nodes->SetLimiter(jPoint,iVar, limiter); - AD::SetPreaccOut(base_nodes->GetLimiter(jPoint)[iVar]); - } - - AD::EndPreacc(); - } - } - - if (LocalMinSolution != NULL) delete [] LocalMinSolution; - if (LocalMaxSolution != NULL) delete [] LocalMaxSolution; - if (GlobalMinSolution != NULL) delete [] GlobalMinSolution; - if (GlobalMaxSolution != NULL) delete [] GlobalMaxSolution; - - } - - /*--- Sharp edges limiter ---*/ - - if (config->GetKind_SlopeLimit() == SHARP_EDGES) { - - /*-- Get limiter parameters from the configuration file ---*/ - - dave = config->GetRefElemLength(); - LimK = config->GetVenkat_LimiterCoeff(); - eps1 = LimK*dave; - eps2 = eps1*eps1*eps1; - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Gradient_i = base_nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = base_nodes->GetGradient_Reconstruction(jPoint); - Coord_i = geometry->node[iPoint]->GetCoord(); - Coord_j = geometry->node[jPoint]->GetCoord(); - - for (iVar = 0; iVar < nVar; iVar++) { - - /*--- Calculate the interface left gradient, delta- (dm) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_j[iDim]-Coord_i[iDim])*Gradient_i[iVar][iDim]; - - /*--- Calculate the interface right gradient, delta+ (dp) ---*/ - - if ( dm > 0.0 ) dp = base_nodes->GetSolution_Max(iPoint, iVar); - else dp = base_nodes->GetSolution_Min(iPoint, iVar); - - /*--- Compute the distance to a sharp edge ---*/ - - SharpEdge_Distance = (geometry->node[iPoint]->GetSharpEdge_Distance() - config->GetAdjSharp_LimiterCoeff()*eps1); - ds = 0.0; - if (SharpEdge_Distance < -eps1) ds = 0.0; - if (fabs(SharpEdge_Distance) <= eps1) ds = 0.5*(1.0+(SharpEdge_Distance/eps1)+(1.0/PI_NUMBER)*sin(PI_NUMBER*SharpEdge_Distance/eps1)); - if (SharpEdge_Distance > eps1) ds = 1.0; - - limiter = ds * ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < base_nodes->GetLimiter(iPoint, iVar)) - base_nodes->SetLimiter(iPoint, iVar, limiter); - - /*-- Repeat for point j on the edge ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_i[iDim]-Coord_j[iDim])*Gradient_j[iVar][iDim]; - - if ( dm > 0.0 ) dp = base_nodes->GetSolution_Max(jPoint,iVar); - else dp = base_nodes->GetSolution_Min(jPoint,iVar); - - /*--- Compute the distance to a sharp edge ---*/ - - SharpEdge_Distance = (geometry->node[jPoint]->GetSharpEdge_Distance() - config->GetAdjSharp_LimiterCoeff()*eps1); - ds = 0.0; - if (SharpEdge_Distance < -eps1) ds = 0.0; - if (fabs(SharpEdge_Distance) <= eps1) ds = 0.5*(1.0+(SharpEdge_Distance/eps1)+(1.0/PI_NUMBER)*sin(PI_NUMBER*SharpEdge_Distance/eps1)); - if (SharpEdge_Distance > eps1) ds = 1.0; - - limiter = ds * ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < base_nodes->GetLimiter(jPoint,iVar)) - base_nodes->SetLimiter(jPoint,iVar, limiter); - - } - } - } - - /*--- Sharp edges limiter ---*/ - - if (config->GetKind_SlopeLimit() == WALL_DISTANCE) { - - /*-- Get limiter parameters from the configuration file ---*/ - - dave = config->GetRefElemLength(); - LimK = config->GetVenkat_LimiterCoeff(); - eps1 = LimK*dave; - eps2 = eps1*eps1*eps1; - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Gradient_i = base_nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = base_nodes->GetGradient_Reconstruction(jPoint); - Coord_i = geometry->node[iPoint]->GetCoord(); - Coord_j = geometry->node[jPoint]->GetCoord(); - - for (iVar = 0; iVar < nVar; iVar++) { - - /*--- Calculate the interface left gradient, delta- (dm) ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_j[iDim]-Coord_i[iDim])*Gradient_i[iVar][iDim]; - - /*--- Calculate the interface right gradient, delta+ (dp) ---*/ - - if ( dm > 0.0 ) dp = base_nodes->GetSolution_Max(iPoint, iVar); - else dp = base_nodes->GetSolution_Min(iPoint, iVar); - - /*--- Compute the distance to a sharp edge ---*/ - - SharpEdge_Distance = (geometry->node[iPoint]->GetWall_Distance() - config->GetAdjSharp_LimiterCoeff()*eps1); - ds = 0.0; - if (SharpEdge_Distance < -eps1) ds = 0.0; - if (fabs(SharpEdge_Distance) <= eps1) ds = 0.5*(1.0+(SharpEdge_Distance/eps1)+(1.0/PI_NUMBER)*sin(PI_NUMBER*SharpEdge_Distance/eps1)); - if (SharpEdge_Distance > eps1) ds = 1.0; - - limiter = ds * ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < base_nodes->GetLimiter(iPoint, iVar)) - base_nodes->SetLimiter(iPoint, iVar, limiter); - - /*-- Repeat for point j on the edge ---*/ - - dm = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - dm += 0.5*(Coord_i[iDim]-Coord_j[iDim])*Gradient_j[iVar][iDim]; - - if ( dm > 0.0 ) dp = base_nodes->GetSolution_Max(jPoint,iVar); - else dp = base_nodes->GetSolution_Min(jPoint,iVar); - - /*--- Compute the distance to a sharp edge ---*/ - - SharpEdge_Distance = (geometry->node[jPoint]->GetWall_Distance() - config->GetAdjSharp_LimiterCoeff()*eps1); - ds = 0.0; - if (SharpEdge_Distance < -eps1) ds = 0.0; - if (fabs(SharpEdge_Distance) <= eps1) ds = 0.5*(1.0+(SharpEdge_Distance/eps1)+(1.0/PI_NUMBER)*sin(PI_NUMBER*SharpEdge_Distance/eps1)); - if (SharpEdge_Distance > eps1) ds = 1.0; - - limiter = ds * ( dp*dp + 2.0*dp*dm + eps2 )/( dp*dp + dp*dm + 2.0*dm*dm + eps2); - - if (limiter < base_nodes->GetLimiter(jPoint,iVar)) - base_nodes->SetLimiter(jPoint,iVar, limiter); - - } - } - } - - /*--- Correct the limiter values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_SOL_2); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LIM_SOL_2); - } - - /*--- Limiter MPI ---*/ - - InitiateComms(geometry, config, SOLUTION_LIMITER); - CompleteComms(geometry, config, SOLUTION_LIMITER); - -#ifdef CODI_REVERSE_TYPE - if (TapeActive) AD::StartRecording(); -#endif + computeLimiters(kindLimiter, this, SOLUTION_LIMITER, PERIODIC_LIM_SOL_1, PERIODIC_LIM_SOL_2, + *geometry, *config, 0, nVar, solution, gradient, solMin, solMax, limiter); } void CSolver::Gauss_Elimination(su2double** A, su2double* rhs, unsigned short nVar) { diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index 45ea2b577193..8d1805a4284c 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -145,10 +145,6 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2 } -void CEulerVariable::SetGradient_PrimitiveZero() { - Gradient_Primitive.storage.setConstant(0.0); -} - bool CEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { bool RightVol = true; diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index c411ba68320b..f809526df41c 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -130,10 +130,6 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci } -void CIncEulerVariable::SetGradient_PrimitiveZero() { - Gradient_Primitive.storage.setConstant(0.0); -} - bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { unsigned long iVar; diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index a5179a0acc1e..acb436bea3fe 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -93,12 +93,6 @@ void CVariable::Restore_BGSSolution_k() { Solution = Solution_BGS_k; } void CVariable::SetResidualSumZero() { Residual_Sum.setConstant(0.0); } -void CVariable::SetAuxVarGradientZero() { Grad_AuxVar.setConstant(0.0); } - -void CVariable::SetGradientZero() { Gradient.storage.setConstant(0.0); } - -void CVariable::SetRmatrixZero() { Rmatrix.storage.setConstant(0.0); } - void CVariable::SetUnd_LaplZero() { Undivided_Laplacian.setConstant(0.0); } void CVariable::SetExternalZero() { External.setConstant(0.0); } diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index a1e0442738f7..a2aec73989fd 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -58,7 +58,7 @@ def main(): naca0012.cfg_dir = "euler/naca0012" naca0012.cfg_file = "inv_NACA0012_Roe.cfg" naca0012.test_iter = 20 - naca0012.test_vals = [-4.056135, -3.565160, 0.336703, 0.021541] #last 4 columns + naca0012.test_vals = [-4.055696, -3.564675, 0.336752, 0.021541] #last 4 columns naca0012.su2_exec = "parallel_computation.py -f" naca0012.timeout = 1600 naca0012.tol = 0.00001 @@ -140,7 +140,7 @@ def main(): cylinder.cfg_dir = "navierstokes/cylinder" cylinder.cfg_file = "lam_cylinder.cfg" cylinder.test_iter = 25 - cylinder.test_vals = [-6.759137, -1.291223, 0.107133, 0.853339] #last 4 columns + cylinder.test_vals = [-6.759136, -1.291221, 0.107150, 0.853257] #last 4 columns cylinder.su2_exec = "parallel_computation.py -f" cylinder.timeout = 1600 cylinder.tol = 0.00001 @@ -173,7 +173,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.493460, -7.672043, -0.000000, 2.085796] #last 4 columns + poiseuille_profile.test_vals = [-12.493462, -7.671815, -0.000000, 2.085796] #last 4 columns poiseuille_profile.su2_exec = "parallel_computation.py -f" poiseuille_profile.timeout = 1600 poiseuille_profile.tol = 0.00001 @@ -232,7 +232,7 @@ def main(): turb_oneram6.cfg_dir = "rans/oneram6" turb_oneram6.cfg_file = "turb_ONERAM6.cfg" turb_oneram6.test_iter = 10 - turb_oneram6.test_vals = [-2.372511, -6.579340, 0.229864, 0.147639] #last 4 columns + turb_oneram6.test_vals = [-2.372346, -6.579371, 0.229867, 0.147637] #last 4 columns turb_oneram6.su2_exec = "parallel_computation.py -f" turb_oneram6.timeout = 3200 turb_oneram6.tol = 0.00001 @@ -243,7 +243,7 @@ def main(): 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.078407, -16.147828, 1.064326, 0.019770] #last 4 columns + turb_naca0012_sa.test_vals = [-12.078482, -16.147828, 1.064326, 0.019770] #last 4 columns turb_naca0012_sa.su2_exec = "parallel_computation.py -f" turb_naca0012_sa.timeout = 3200 turb_naca0012_sa.tol = 0.00001 @@ -509,7 +509,7 @@ def main(): contadj_naca0012.cfg_dir = "cont_adj_euler/naca0012" contadj_naca0012.cfg_file = "inv_NACA0012.cfg" contadj_naca0012.test_iter = 5 - contadj_naca0012.test_vals = [-9.783199, -15.190764, 3.0092e-01, 1.9552e-02] #last 4 columns + contadj_naca0012.test_vals = [-9.300829, -14.583240, 0.300920, 0.019552] #last 4 columns contadj_naca0012.su2_exec = "parallel_computation.py -f" contadj_naca0012.timeout = 1600 contadj_naca0012.tol = 0.00001 @@ -631,7 +631,7 @@ def main(): 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.906256, 1.337943, 6.052217, 2.395627] #last 4 columns + turb_naca0012_1c.test_vals = [-4.906243, 1.337959, 6.052221, 2.395629] #last 4 columns turb_naca0012_1c.su2_exec = "parallel_computation.py -f" turb_naca0012_1c.timeout = 1600 turb_naca0012_1c.tol = 0.00001 @@ -744,7 +744,7 @@ def main(): square_cylinder.cfg_dir = "unsteady/square_cylinder" square_cylinder.cfg_file = "turb_square.cfg" square_cylinder.test_iter = 3 - square_cylinder.test_vals = [-1.162621, 0.066395, 1.399789, 2.220408] #last 4 columns + square_cylinder.test_vals = [-1.162660, 0.066413, 1.399789, 2.220408] #last 4 columns square_cylinder.su2_exec = "parallel_computation.py -f" square_cylinder.timeout = 1600 square_cylinder.tol = 0.00001 @@ -870,7 +870,7 @@ def main(): 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.574121, 5.820564, 96.994080, 0.062865] #last 4 columns + transonic_stator.test_vals = [-0.574097, 5.820577, 96.994270, 0.062868] #last 4 columns transonic_stator.su2_exec = "parallel_computation.py -f" transonic_stator.timeout = 1600 transonic_stator.new_output = False @@ -882,7 +882,7 @@ def main(): 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 = [-2.125313, 3.014645, 5.300089, 0.003103] #last 4 columns + transonic_stator_rst.test_vals = [-6.797150, -0.755136, 5.007406, 0.0029491] #last 4 columns transonic_stator_rst.su2_exec = "parallel_computation.py -f" transonic_stator_rst.timeout = 1600 transonic_stator_rst.new_output = False @@ -898,7 +898,7 @@ def main(): 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.188747, -10.631530] #last 4 columns + uniform_flow.test_vals = [5.000000, 0.000000, -0.188747, -10.631534] #last 4 columns uniform_flow.su2_exec = "parallel_computation.py -f" uniform_flow.timeout = 1600 uniform_flow.tol = 0.000001 @@ -911,7 +911,7 @@ def main(): 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.398126, 0.353070, 0.405696] #last 4 columns + channel_2D.test_vals = [2.000000, 0.000000, 0.398157, 0.353079, 0.405679] #last 4 columns channel_2D.su2_exec = "parallel_computation.py -f" channel_2D.timeout = 100 channel_2D.tol = 0.00001 @@ -924,7 +924,7 @@ def main(): 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.620109, 0.505162, 0.415445] #last 4 columns + channel_3D.test_vals = [2.000000, 0.000000, 0.620137, 0.505221, 0.415605] #last 4 columns channel_3D.su2_exec = "parallel_computation.py -f" channel_3D.timeout = 1600 channel_3D.tol = 0.00001 @@ -963,7 +963,7 @@ def main(): supersonic_vortex_shedding.cfg_dir = "sliding_interface/supersonic_vortex_shedding" supersonic_vortex_shedding.cfg_file = "sup_vor_shed_WA.cfg" supersonic_vortex_shedding.test_iter = 5 - supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.244192, 1.644524] #last 4 columns + supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.228298, 1.648277] #last 4 columns supersonic_vortex_shedding.su2_exec = "parallel_computation.py -f" supersonic_vortex_shedding.timeout = 1600 supersonic_vortex_shedding.tol = 0.00001 @@ -1109,7 +1109,7 @@ def main(): pywrapper_naca0012.cfg_dir = "euler/naca0012" pywrapper_naca0012.cfg_file = "inv_NACA0012_Roe.cfg" pywrapper_naca0012.test_iter = 100 - pywrapper_naca0012.test_vals = [-7.306809, -6.740117, 0.333437, 0.021233] #last 4 columns + pywrapper_naca0012.test_vals = [-7.298733, -6.733533, 0.333456, 0.021233] #last 4 columns pywrapper_naca0012.su2_exec = "mpirun -np 2 SU2_CFD.py --parallel -f" pywrapper_naca0012.timeout = 1600 pywrapper_naca0012.tol = 0.00001 @@ -1131,7 +1131,7 @@ def main(): pywrapper_square_cylinder.cfg_dir = "unsteady/square_cylinder" pywrapper_square_cylinder.cfg_file = "turb_square.cfg" pywrapper_square_cylinder.test_iter = 3 - pywrapper_square_cylinder.test_vals = [-1.162621, 0.066395, 1.399789, 2.220408] #last 4 columns + pywrapper_square_cylinder.test_vals = [-1.162660, 0.066413, 1.399789, 2.220408] #last 4 columns pywrapper_square_cylinder.su2_exec = "mpirun -np 2 SU2_CFD.py --parallel -f" pywrapper_square_cylinder.timeout = 1600 pywrapper_square_cylinder.tol = 0.00001 diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index c73e215f831d..68b7ebce3167 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -166,7 +166,7 @@ def main(): discadj_cylinder.cfg_dir = "disc_adj_rans/cylinder" discadj_cylinder.cfg_file = "cylinder.cfg" discadj_cylinder.test_iter = 9 - discadj_cylinder.test_vals = [3.746900, -1.544893, -8.3447e-03, 1.3808e-05] #last 4 columns + discadj_cylinder.test_vals = [3.746909, -1.544883, -0.008321, 0.000014] #last 4 columns discadj_cylinder.su2_exec = "parallel_computation.py -f" discadj_cylinder.timeout = 1600 discadj_cylinder.tol = 0.00001 @@ -182,7 +182,7 @@ def main(): discadj_cylinder.cfg_dir = "disc_adj_rans/cylinder" discadj_cylinder.cfg_file = "cylinder_Windowing_AD.cfg" discadj_cylinder.test_iter = 9 - discadj_cylinder.test_vals = [3.003855] #last column + discadj_cylinder.test_vals = [3.004406] #last column discadj_cylinder.su2_exec = "parallel_computation.py -f" discadj_cylinder.timeout = 1600 discadj_cylinder.tol = 0.00001 @@ -215,7 +215,7 @@ def main(): discadj_DT_1ST_cylinder.cfg_dir = "disc_adj_rans/cylinder_DT_1ST" discadj_DT_1ST_cylinder.cfg_file = "cylinder.cfg" discadj_DT_1ST_cylinder.test_iter = 9 - discadj_DT_1ST_cylinder.test_vals = [3.698165, -1.607052, -2.2503e-03, 2.7212e-05] #last 4 columns + discadj_DT_1ST_cylinder.test_vals = [3.698168, -1.607050, -0.002159, 0.000028] #last 4 columns discadj_DT_1ST_cylinder.su2_exec = "parallel_computation.py -f" discadj_DT_1ST_cylinder.timeout = 1600 discadj_DT_1ST_cylinder.tol = 0.00001 @@ -247,7 +247,7 @@ def main(): discadj_trans_stator.cfg_dir = "disc_adj_turbomachinery/transonic_stator_2D" discadj_trans_stator.cfg_file = "transonic_stator.cfg" discadj_trans_stator.test_iter = 79 - discadj_trans_stator.test_vals = [79.000000, -1.927296, -1.401205] #last 4 columns + discadj_trans_stator.test_vals = [79.000000, -1.927271, -1.401487] #last 4 columns discadj_trans_stator.su2_exec = "parallel_computation.py -f" discadj_trans_stator.timeout = 1600 discadj_trans_stator.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 3911b8cd474d..332cf3977aa7 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -57,7 +57,7 @@ def main(): naca0012.cfg_dir = "euler/naca0012" naca0012.cfg_file = "inv_NACA0012_Roe.cfg" naca0012.test_iter = 20 - naca0012.test_vals = [-4.024286, -3.515166, 0.339392, 0.022217] #last 4 columns + naca0012.test_vals = [-4.023999, -3.515034, 0.339426, 0.022217] #last 4 columns naca0012.su2_exec = "SU2_CFD" naca0012.timeout = 1600 naca0012.new_output= True @@ -146,7 +146,7 @@ def main(): cylinder.cfg_dir = "navierstokes/cylinder" cylinder.cfg_file = "lam_cylinder.cfg" cylinder.test_iter = 25 - cylinder.test_vals = [-6.765432, -1.297428, 0.019508, 0.310040] #last 4 columns + cylinder.test_vals = [-6.765430, -1.297426, 0.019508, 0.310015] #last 4 columns cylinder.su2_exec = "SU2_CFD" cylinder.new_output = True cylinder.timeout = 1600 @@ -182,7 +182,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.494724, -7.712336, -0.000000, 2.085796] #last 4 columns + poiseuille_profile.test_vals = [-12.494733, -7.712320, -0.000000, 2.085796] #last 4 columns poiseuille_profile.su2_exec = "SU2_CFD" poiseuille_profile.new_output = True poiseuille_profile.timeout = 1600 @@ -245,7 +245,7 @@ def main(): turb_oneram6.cfg_dir = "rans/oneram6" turb_oneram6.cfg_file = "turb_ONERAM6.cfg" turb_oneram6.test_iter = 10 - turb_oneram6.test_vals = [-2.372511, -6.579339, 0.229864, 0.147639]#last 4 columns + turb_oneram6.test_vals = [-2.372347, -6.579371, 0.229867, 0.147638]#last 4 columns turb_oneram6.su2_exec = "SU2_CFD" turb_oneram6.new_output = True turb_oneram6.timeout = 3200 @@ -257,7 +257,7 @@ def main(): 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.075873, -16.146770, 1.064326, 0.019770] #last 4 columns + turb_naca0012_sa.test_vals = [-12.075973, -16.146770, 1.064326, 0.019770] #last 4 columns turb_naca0012_sa.su2_exec = "SU2_CFD" turb_naca0012_sa.new_output = True turb_naca0012_sa.timeout = 3200 @@ -527,7 +527,7 @@ def main(): schubauer_klebanoff_transition.cfg_file = "transitional_BC_model_ConfigFile.cfg" schubauer_klebanoff_transition.test_iter = 10 schubauer_klebanoff_transition.new_output = True - schubauer_klebanoff_transition.test_vals = [-8.029756, -14.268351, 0.000053, 0.007986] #last 4 columns + schubauer_klebanoff_transition.test_vals = [-8.029786, -14.268351, 0.000053, 0.007986] #last 4 columns schubauer_klebanoff_transition.su2_exec = "SU2_CFD" schubauer_klebanoff_transition.timeout = 1600 schubauer_klebanoff_transition.tol = 0.00001 @@ -542,7 +542,7 @@ def main(): contadj_naca0012.cfg_dir = "cont_adj_euler/naca0012" contadj_naca0012.cfg_file = "inv_NACA0012.cfg" contadj_naca0012.test_iter = 5 - contadj_naca0012.test_vals = [-9.787554, -15.192510, 3.0092e-01, 1.9552e-02] #last 4 columns + contadj_naca0012.test_vals = [-9.289565, -14.563861, 0.300920, 0.019552] #last 4 columns contadj_naca0012.su2_exec = "SU2_CFD" contadj_naca0012.new_output = True contadj_naca0012.timeout = 1600 @@ -710,7 +710,7 @@ def main(): 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.002978, 1.312367, 6.085205, 2.413463] #last 4 columns + turb_naca0012_p1c1.test_vals = [-5.003083, 1.312255, 6.085205, 2.413463] #last 4 columns turb_naca0012_p1c1.su2_exec = "SU2_CFD" turb_naca0012_p1c1.new_output = True turb_naca0012_p1c1.timeout = 1600 @@ -794,7 +794,7 @@ def main(): square_cylinder.cfg_dir = "unsteady/square_cylinder" square_cylinder.cfg_file = "turb_square.cfg" square_cylinder.test_iter = 3 - square_cylinder.test_vals = [-1.162581, 0.066412, 1.399788, 2.220411] #last 4 columns + square_cylinder.test_vals = [-1.162561, 0.066414, 1.399788, 2.220411] #last 4 columns square_cylinder.su2_exec = "SU2_CFD" square_cylinder.timeout = 1600 square_cylinder.tol = 0.00001 @@ -939,7 +939,7 @@ def main(): 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.560481, 5.828675, 96.593740, 0.062507] #last 4 columns + transonic_stator.test_vals = [-0.560463, 5.828675, 96.594070, 0.062505] #last 4 columns transonic_stator.su2_exec = "SU2_CFD" transonic_stator.new_output = False transonic_stator.timeout = 1600 @@ -951,7 +951,7 @@ def main(): 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 = [-2.127925, 2.936696, 5.298641, 0.003102] #last 4 columns + transonic_stator_rst.test_vals = [-6.781212, -0.757588, 5.007406, 0.0029491] #last 4 columns transonic_stator_rst.su2_exec = "SU2_CFD" transonic_stator_rst.new_output = False transonic_stator_rst.timeout = 1600 @@ -968,7 +968,7 @@ def main(): uniform_flow.cfg_dir = "sliding_interface/uniform_flow" uniform_flow.cfg_file = "uniform_NN.cfg" uniform_flow.test_iter = 2 - uniform_flow.test_vals = [2.000000, 0.000000, -0.205134, -13.254849] #last 4 columns + uniform_flow.test_vals = [2.000000, 0.000000, -0.205134, -13.254129] #last 4 columns uniform_flow.su2_exec = "SU2_CFD" uniform_flow.timeout = 1600 uniform_flow.tol = 0.000001 @@ -981,7 +981,7 @@ def main(): 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.397956, 0.353078, 0.405707] #last 4 columns + channel_2D.test_vals = [2.000000, 0.000000, 0.398070, 0.353082, 0.405729] #last 4 columns channel_2D.su2_exec = "SU2_CFD" channel_2D.timeout = 100 channel_2D.tol = 0.00001 @@ -994,7 +994,7 @@ def main(): channel_3D.cfg_dir = "sliding_interface/channel_3D" channel_3D.cfg_file = "channel_3D_WA.cfg" channel_3D.test_iter = 1 - channel_3D.test_vals = [1.000000, 0.000000, 0.661392, 0.769773, 0.696154] #last 5 columns + channel_3D.test_vals = [1.000000, 0.000000, 0.661403, 0.769795, 0.696074] #last 5 columns channel_3D.su2_exec = "SU2_CFD" channel_3D.timeout = 1600 channel_3D.tol = 0.00001 @@ -1033,7 +1033,7 @@ def main(): supersonic_vortex_shedding.cfg_dir = "sliding_interface/supersonic_vortex_shedding" supersonic_vortex_shedding.cfg_file = "sup_vor_shed_WA.cfg" supersonic_vortex_shedding.test_iter = 5 - supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.240131, 1.645383] #last 4 columns + supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.228735, 1.648284] #last 4 columns supersonic_vortex_shedding.su2_exec = "SU2_CFD" supersonic_vortex_shedding.timeout = 1600 supersonic_vortex_shedding.tol = 0.00001 @@ -1562,7 +1562,7 @@ def main(): pywrapper_naca0012.cfg_dir = "euler/naca0012" pywrapper_naca0012.cfg_file = "inv_NACA0012_Roe.cfg" pywrapper_naca0012.test_iter = 20 - pywrapper_naca0012.test_vals = [-4.024286, -3.515166, 0.339392, 0.022217] #last 4 columns + pywrapper_naca0012.test_vals = [-4.023999, -3.515034, 0.339426, 0.022217] #last 4 columns pywrapper_naca0012.su2_exec = "SU2_CFD.py -f" pywrapper_naca0012.new_output = True pywrapper_naca0012.timeout = 1600 @@ -1588,7 +1588,7 @@ def main(): pywrapper_square_cylinder.cfg_dir = "unsteady/square_cylinder" pywrapper_square_cylinder.cfg_file = "turb_square.cfg" pywrapper_square_cylinder.test_iter = 3 - pywrapper_square_cylinder.test_vals = [-1.162581, 0.066412, 1.399788, 2.220411] #last 4 columns + pywrapper_square_cylinder.test_vals = [-1.162561, 0.066414, 1.399788, 2.220411] #last 4 columns pywrapper_square_cylinder.su2_exec = "SU2_CFD.py -f" pywrapper_square_cylinder.timeout = 1600 pywrapper_square_cylinder.tol = 0.00001 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index ab2d9f613812..841249a6da3c 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -166,7 +166,7 @@ def main(): discadj_cylinder.cfg_dir = "disc_adj_rans/cylinder" discadj_cylinder.cfg_file = "cylinder.cfg" discadj_cylinder.test_iter = 9 - discadj_cylinder.test_vals = [3.746904, -1.544886, -0.008345, 0.000014] #last 4 columns + discadj_cylinder.test_vals = [3.746909, -1.544883, -0.008321, 0.000014] #last 4 columns discadj_cylinder.su2_exec = "SU2_CFD_AD" discadj_cylinder.timeout = 1600 discadj_cylinder.tol = 0.00001 @@ -182,7 +182,7 @@ def main(): discadj_DT_1ST_cylinder.cfg_dir = "disc_adj_rans/cylinder_DT_1ST" discadj_DT_1ST_cylinder.cfg_file = "cylinder.cfg" discadj_DT_1ST_cylinder.test_iter = 9 - discadj_DT_1ST_cylinder.test_vals = [3.698165, -1.607052, -2.2500e-03, 2.7211e-05] #last 4 columns + discadj_DT_1ST_cylinder.test_vals = [3.698168, -1.607050, -0.002159, 0.000028] #last 4 columns discadj_DT_1ST_cylinder.su2_exec = "SU2_CFD_AD" discadj_DT_1ST_cylinder.timeout = 1600 discadj_DT_1ST_cylinder.tol = 0.00001 diff --git a/TestCases/turbomachinery/transonic_stator_2D/transonic_stator_rst.cfg b/TestCases/turbomachinery/transonic_stator_2D/transonic_stator_rst.cfg index d8d209589124..404834b6cf9e 100644 --- a/TestCases/turbomachinery/transonic_stator_2D/transonic_stator_rst.cfg +++ b/TestCases/turbomachinery/transonic_stator_2D/transonic_stator_rst.cfg @@ -243,7 +243,7 @@ CONV_NUM_METHOD_FLOW= ROE MUSCL_FLOW= YES % % Slope limiter (VENKATAKRISHNAN, VAN_ALBADA_EDGE) -SLOPE_LIMITER_FLOW= VENKATAKRISHNAN +SLOPE_LIMITER_FLOW= VENKATAKRISHNAN_WANG % % Entropy fix coefficient (0.0 implies no entropy fixing, 1.0 implies scalar artificial dissipation, 0.001 default) ENTROPY_FIX_COEFF= 0.03