diff --git a/Common/include/containers/container_decorators.hpp b/Common/include/containers/container_decorators.hpp index a9e66857eb3..331e87b9a0f 100644 --- a/Common/include/containers/container_decorators.hpp +++ b/Common/include/containers/container_decorators.hpp @@ -62,6 +62,33 @@ class CMatrixView { const Scalar* operator[](Index i) const noexcept { return &m_ptr[i * m_cols]; } const Scalar& operator()(Index i, Index j) const noexcept { return m_ptr[i * m_cols + j]; } + /*! + * \brief Return copy of data in a static size container (see C2DContainer::get). + * \param[in] i - Row of the view (e.g. point index, whole-mesh usage). + * \param[in] start - Starting column to copy the data (amount determined by container size). + */ + template + StaticContainer get(Index i, Index start = 0) const noexcept { + constexpr size_t Size = StaticContainer::StaticSize; + static_assert(Size, "This method requires a static output type."); + StaticContainer ret; + for (size_t k = 0; k < Size; ++k) ret.data()[k] = m_ptr[i * m_cols + start + k]; + return ret; + } + + /*! + * \brief SIMD gather version of get, one row per lane. + */ + template + StaticContainer get(simd::Array i, Index start = 0) const noexcept { + constexpr size_t Size = StaticContainer::StaticSize; + static_assert(Size, "This method requires a static output type."); + StaticContainer ret; + for (size_t lane = 0; lane < N; ++lane) + for (size_t k = 0; k < Size; ++k) ret.data()[k][lane] = m_ptr[i[lane] * m_cols + start + k]; + return ret; + } + template ::value> = 0> Scalar* operator[](Index i) noexcept { return &m_ptr[i * m_cols]; diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 384b0bf6212..85f3b152be3 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -1055,6 +1055,80 @@ class CSysMatrix { SetBlocks(iEdge, block_i, block_j, -1); } + /*! + * \brief Set the four blocks of an edge, for fluxes whose i and j contributions are independent. + * \note The diagonal blocks are accumulated, the off-diagonal blocks are set. + */ + template + inline void SetBlocks(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, const MatrixType& jac_ii, + const MatrixType& jac_ij, const MatrixType& jac_ji, const MatrixType& jac_jj, + OtherType mask = 1) { + const auto blkSz = nVar * nEqn; + auto* bii = &mat.d[iPoint * blkSz]; + auto* bjj = &mat.d[jPoint * blkSz]; + unsigned long iVar, jVar, offset = 0; + + if (quantized_mode) { + ScalarType bij_buf[MAXNVAR * MAXNVAR], bji_buf[MAXNVAR * MAXNVAR]; + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++, ++offset) { + bii[offset] += PassiveAssign(jac_ii[iVar][jVar] * mask); + bjj[offset] += PassiveAssign(jac_jj[iVar][jVar] * mask); + bij_buf[offset] = PassiveAssign(jac_ij[iVar][jVar] * mask); + bji_buf[offset] = PassiveAssign(jac_ji[iVar][jVar] * mask); + } + QuantizeBlock(bij_buf, &q_scale.u[iEdge * nVar], &q_blocks.u[iEdge * blkSz]); + const auto k_l = edge_ptr_l[iEdge]; + QuantizeBlock(bji_buf, &q_scale.l[k_l * nVar], &q_blocks.l[k_l * blkSz]); + return; + } + + auto* bij = &mat.u[iEdge * blkSz]; + auto* bji = &mat.l[edge_ptr_l[iEdge] * blkSz]; + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nEqn; jVar++) { + bii[offset] += PassiveAssign(jac_ii[iVar][jVar] * mask); + bjj[offset] += PassiveAssign(jac_jj[iVar][jVar] * mask); + bij[offset] = PassiveAssign(jac_ij[iVar][jVar] * mask); + bji[offset] = PassiveAssign(jac_ji[iVar][jVar] * mask); + ++offset; + } + } + } + + /*! + * \brief Set the off-diagonal blocks of an edge, the diagonal being assembled elsewhere. + */ + template + inline void SetOffDiagBlocks(unsigned long iEdge, const MatrixType& jac_ij, const MatrixType& jac_ji, + OtherType mask = 1) { + const auto blkSz = nVar * nEqn; + unsigned long iVar, jVar, offset = 0; + + if (quantized_mode) { + ScalarType bij_buf[MAXNVAR * MAXNVAR], bji_buf[MAXNVAR * MAXNVAR]; + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nEqn; jVar++, ++offset) { + bij_buf[offset] = PassiveAssign(jac_ij[iVar][jVar] * mask); + bji_buf[offset] = PassiveAssign(jac_ji[iVar][jVar] * mask); + } + QuantizeBlock(bij_buf, &q_scale.u[iEdge * nVar], &q_blocks.u[iEdge * blkSz]); + const auto k_l = edge_ptr_l[iEdge]; + QuantizeBlock(bji_buf, &q_scale.l[k_l * nVar], &q_blocks.l[k_l * blkSz]); + return; + } + + auto* bij = &mat.u[iEdge * blkSz]; + auto* bji = &mat.l[edge_ptr_l[iEdge] * blkSz]; + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nEqn; jVar++) { + bij[offset] = PassiveAssign(jac_ij[iVar][jVar] * mask); + bji[offset] = PassiveAssign(jac_ji[iVar][jVar] * mask); + ++offset; + } + } + } + /*! * \brief SIMD version, does the update for multiple edges. * \note Nothing is updated if the mask is 0. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index c21116c6f3e..5851b10a0b4 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4352,6 +4352,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } /* --- Check for NEMO compatibility issues ---*/ + if (nemo && Kind_Turb_Model != TURB_MODEL::NONE) { + SU2_MPI::Error("A turbulence model is not yet available for the NEMO solver.", CURRENT_FUNCTION); + } + if (Kind_FluidModel == SU2_NONEQ && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::SUTHERLAND && Kind_TransCoeffModel != TRANSCOEFFMODEL::GUPTAYOS) ) { SU2_MPI::Error("Transport model not available for NEMO solver using SU2TCLIB. Please use the WILKE, SUTHERLAND or GUPTAYOS transport model instead.", CURRENT_FUNCTION); } diff --git a/SU2_CFD/include/numerics/heat.hpp b/SU2_CFD/include/numerics/heat.hpp deleted file mode 100644 index 1e9cbf81505..00000000000 --- a/SU2_CFD/include/numerics/heat.hpp +++ /dev/null @@ -1,107 +0,0 @@ -/*! - * \file heat.hpp - * \brief Declarations of numerics classes for heat transfer problems. - * \author F. Palacios, T. Economon - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, 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 . - */ - -#pragma once - -#include "scalar/scalar_diffusion.hpp" -#include "scalar/scalar_convection.hpp" -#include "../variables/CIncEulerVariable.hpp" - -/*! - * \class CUpwSca_Heat - * \brief Class for doing a scalar upwind solver for the heat convection equation. - * \ingroup ConvDiscr - * \author O. Burghardt. - * \version 8.5.0 "Harrier" - */ -class CUpwSca_Heat final : public CUpwScalar> { - public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] config - Definition of the particular problem. - */ - CUpwSca_Heat(unsigned short val_nDim, const CConfig *config) - : CUpwScalar>(val_nDim, 1, config) {} - - private: - /*! - * \brief Adds extra variables to AD - */ - void ExtraADPreaccIn(void) override {} - - /*! - * \brief Heat-specific specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - Flux[0] = a0 * ScalarVar_i[0] + a1 * ScalarVar_j[0]; - Jacobian_i[0][0] = a0; - Jacobian_j[0][0] = a1; - } -}; - -/*! - * \class CAvgGrad_Heat - * \brief Class for computing viscous term using average of gradients without correction (heat equation). - * \ingroup ViscDiscr - * \author O. Burghardt. - * \version 8.5.0 "Harrier" - */ -class CAvgGrad_Heat final : public CAvgGrad_Scalar { - public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] correct - Whether to correct the gradient. - */ - CAvgGrad_Heat(unsigned short val_nDim, const CConfig *config, bool correct) - : CAvgGrad_Scalar(val_nDim, 1, correct, config) {} - - private: - /*! - * \brief Adds extra variables to AD - */ - void ExtraADPreaccIn(void) override { - AD::SetPreaccIn(*Diffusion_Coeff_i, *Diffusion_Coeff_j); - } - - /*! - * \brief Heat-specific specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - const su2double Thermal_Diffusivity_Mean = 0.5 * (*Diffusion_Coeff_i + *Diffusion_Coeff_j); - - Flux[0] = Thermal_Diffusivity_Mean * Proj_Mean_GradScalarVar[0]; - - /*--- Use TSL for Jacobians. ---*/ - Jacobian_i[0][0] = -Thermal_Diffusivity_Mean * proj_vector_ij; - Jacobian_j[0][0] = Thermal_Diffusivity_Mean * proj_vector_ij; - } -}; diff --git a/SU2_CFD/include/numerics/heat_edge_flux.hpp b/SU2_CFD/include/numerics/heat_edge_flux.hpp new file mode 100644 index 00000000000..213d168d97f --- /dev/null +++ b/SU2_CFD/include/numerics/heat_edge_flux.hpp @@ -0,0 +1,89 @@ +/*! + * \file heat_edge_flux.hpp + * \brief Heat transport as a third-layer scalar flux, see numerics/scalar/scalar_edge_flux.hpp. + * \author P. Gomes + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, 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 . + */ + +#pragma once + +#include "scalar/scalar_edge_flux.hpp" + +/*! + * \class CScalarFlux_Heat + * \ingroup ViscDiscr + * \brief Convection and diffusion of temperature, non-conservative, with a diagonal + * (single-equation) diffusion coefficient. + * \note The temperature has no notion of density weighting, so Conservative is false and the + * inherited CUpwScalarFlux::finalizeFlux, flux(0) = a0*phi_i(0) + a1*phi_j(0), is exactly + * the model's whole convective term; no override is needed. + * \note The solver runs in two modes, weakly-coupled energy equation on a fluid zone or standalone + * conduction on a solid one (CHeatSolver::flow); the diffusion coefficient is the flow's + * thermal conductivity over specific heat, plus a turbulent contribution, in the former, and + * the configured constant thermal diffusivity in the latter. EdgeSide::flowNodes is null in + * the solid case, so this is the only place that may read it, and only when flow is set. + */ +template +class CScalarFlux_Heat final + : public CUpwScalarBase, FlowIndices, nDim, nVar> { + public: + static constexpr bool Conservative = false; + static constexpr bool DiagonalDiffusion = true; + + using Base = CUpwScalarBase; + using Int = typename Base::Int; + + explicit CScalarFlux_Heat(const CConfig& config) + : Base(config), + flow(config.GetFluidProblem()), + prandtlTurb(config.GetPrandtl_Turb()), + constDiffusivity(config.GetThermalDiffusivity()) {} + + /*! + * \brief Thermal diffusivity, an i/j average, identical for both edge sides (TSL Jacobian). + */ + template + FORCEINLINE CPair> coefficients(const FlowIndices& idx, Int iPoint, + const EdgeSide& side_i, Int jPoint, + const EdgeSide& side_j, + const CPair&) const { + Vector D; + if (flow) { + const Double k_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.ThermalConductivity()); + const Double cp_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.CpTotal()); + const Double muT_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.EddyViscosity()); + const Double k_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.ThermalConductivity()); + const Double cp_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.CpTotal()); + const Double muT_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.EddyViscosity()); + D(0) = 0.5 * (k_i / cp_i + muT_i / prandtlTurb + k_j / cp_j + muT_j / prandtlTurb); + } else { + D(0) = constDiffusivity; + } + return {D, D}; + } + + private: + const bool flow; + const su2double prandtlTurb; + const su2double constDiffusivity; +}; diff --git a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp deleted file mode 100644 index e40749a5c1d..00000000000 --- a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp +++ /dev/null @@ -1,150 +0,0 @@ -/*! - * \file scalar_convection.hpp - * \brief Declarations of numerics classes for discretization of - * convective fluxes in scalar problems. - * \author F. Palacios, T. Economon - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, 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 . - */ - -#pragma once - -#include "../CNumerics.hpp" - -/*! - * \class CUpwScalar - * \brief Template class for scalar upwind fluxes between nodes i and j. - * \details This class serves as a template for the scalar upwinding residual - * classes. The general structure of a scalar upwinding calculation is the - * same for many different models, which leads to a lot of repeated code. - * By using the template design pattern, these sections of repeated code are - * moved to this shared base class, and the specifics of each model - * are implemented by derived classes. In order to add a new residual - * calculation for a convection residual, extend this class and implement - * the pure virtual functions with model-specific behavior. - * \ingroup ConvDiscr - * \author C. Pederson, A. Bueno., and A. Campos. - */ -template -class CUpwScalar : public CNumerics { - protected: - enum : unsigned short {MAXNVAR = 8}; - - const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ - su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0. */ - su2double a1 = 0.0; /*!< \brief The minimum of the face-normal velocity and 0. */ - su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */ - su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */ - su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */ - su2double JacobianBuffer[2*MAXNVAR*MAXNVAR]; /*!< \brief Static storage for the two Jacobians. */ - - const bool incompressible = false, dynamic_grid = false; - - /*! - * \brief A pure virtual function. Derived classes must use it to register the additional - * variables they use as preaccumulation inputs, e.g. the density for SST. - */ - virtual void ExtraADPreaccIn() = 0; - - /*! - * \brief Model-specific steps in the ComputeResidual method, derived classes - * compute the Flux and its Jacobians via this method. - * \param[in] config - Definition of the particular problem. - */ - virtual void FinishResidualCalc(const CConfig* config) = 0; - - public: - /*! - * \brief Constructor of the class. - * \param[in] ndim - Number of dimensions of the problem. - * \param[in] nvar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CUpwScalar(unsigned short ndim, unsigned short nvar, const CConfig* config) - : CNumerics(ndim, nvar, config), - idx(ndim, config->GetnSpecies()), - incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE), - dynamic_grid(config->GetDynamic_Grid()) { - if (nVar > MAXNVAR) { - SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION); - } - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = &JacobianBuffer[iVar * nVar]; - Jacobian_j[iVar] = &JacobianBuffer[iVar * nVar + MAXNVAR * MAXNVAR]; - } - - /*--- Initialize the JacobianBuffer to zero. ---*/ - for (unsigned short iVar = 0; iVar < 2*MAXNVAR*MAXNVAR; iVar++) { - JacobianBuffer[iVar] = 0.0; - } - } - - /*! - * \brief Compute the scalar upwind flux between two nodes i and j. - * \param[in] config - Definition of the particular problem. - * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. - */ - CNumerics::ResidualType<> ComputeResidual(const CConfig* config) final { - AD::StartPreacc(); - AD::SetPreaccIn(Normal, nDim); - AD::SetPreaccIn(ScalarVar_i, nVar); - AD::SetPreaccIn(ScalarVar_j, nVar); - if (dynamic_grid) { - AD::SetPreaccIn(GridVel_i, nDim); - AD::SetPreaccIn(GridVel_j, nDim); - } - AD::SetPreaccIn(&V_i[idx.Velocity()], nDim); - AD::SetPreaccIn(&V_j[idx.Velocity()], nDim); - AD::SetPreaccIn(V_i[idx.Density()]); - AD::SetPreaccIn(V_j[idx.Density()]); - AD::SetPreaccIn(MassFlux); - - ExtraADPreaccIn(); - - if (bounded_scalar) { - a0 = fmax(0.0, MassFlux) / V_i[idx.Density()]; - a1 = fmin(0.0, MassFlux) / V_j[idx.Density()]; - } else { - su2double q_ij = 0.0; - if (dynamic_grid) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) { - su2double Velocity_i = V_i[iDim + idx.Velocity()] - GridVel_i[iDim]; - su2double Velocity_j = V_j[iDim + idx.Velocity()] - GridVel_j[iDim]; - q_ij += 0.5 * (Velocity_i + Velocity_j) * Normal[iDim]; - } - } else { - for (unsigned short iDim = 0; iDim < nDim; iDim++) { - q_ij += 0.5 * (V_i[iDim + idx.Velocity()] + V_j[iDim + idx.Velocity()]) * Normal[iDim]; - } - } - a0 = fmax(0.0, q_ij); - a1 = fmin(0.0, q_ij); - } - - FinishResidualCalc(config); - - AD::SetPreaccOut(Flux, nVar); - AD::EndPreacc(); - - return ResidualType<>(Flux, Jacobian_i, Jacobian_j); - } -}; diff --git a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp deleted file mode 100644 index 580c0bedcef..00000000000 --- a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp +++ /dev/null @@ -1,155 +0,0 @@ -/*! - * \file scalar_diffusion.hpp - * \brief Declarations of numerics classes for discretization of - * viscous fluxes in scalar problems. - * \author F. Palacios, T. Economon - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, 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 . - */ - -#pragma once - -#include "../CNumerics.hpp" - -/*! - * \class CNoFlowIndices - * \brief Dummy flow indices class to use CAvgGrad_Scalar when flow variables are not available. - * For example, solid heat transfer problems. - */ -struct CNoFlowIndices { - CNoFlowIndices(int, int) {} - inline int Density() const { return 0; } - inline int LaminarViscosity() const { return 0; } - inline int EddyViscosity() const { return 0; } -}; - -/*! - * \class CAvgGrad_Scalar - * \brief Template class for computing viscous residual of scalar values - * \details This class serves as a template for the scalar viscous residual - * classes. The general structure of a viscous residual calculation is the - * same for many different models, which leads to a lot of repeated code. - * By using the template design pattern, these sections of repeated code are - * moved to a shared base class, and the specifics of each model - * are implemented by derived classes. In order to add a new residual - * calculation for a viscous residual, extend this class and implement - * the pure virtual functions with model-specific behavior. - * \ingroup ViscDiscr - * \author C. Pederson, A. Bueno, and F. Palacios - */ -template -class CAvgGrad_Scalar : public CNumerics { - protected: - enum : unsigned short {MAXNVAR = 8}; - - const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ - su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */ - su2double proj_vector_ij = 0.0; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ - su2double Flux[MAXNVAR] = {0.0}; /*!< \brief Final result, diffusive flux/residual. */ - su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */ - su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */ - su2double JacobianBuffer[2*MAXNVAR*MAXNVAR];/*!< \brief Static storage for the two Jacobians. */ - - const bool correct_gradient = false, incompressible = false; - - /*! - * \brief A pure virtual function; Adds any extra variables to AD - */ - virtual void ExtraADPreaccIn() = 0; - - /*! - * \brief Model-specific steps in the ComputeResidual method, derived classes - * should compute the Flux and Jacobians (i/j) inside this method. - * \param[in] config - Definition of the particular problem. - */ - virtual void FinishResidualCalc(const CConfig* config) = 0; - - public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] correct_gradient - Whether to correct gradient for skewness. - * \param[in] config - Definition of the particular problem. - */ - CAvgGrad_Scalar(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, - const CConfig* config) - : CNumerics(val_nDim, val_nVar, config), - idx(val_nDim, config->GetnSpecies()), - correct_gradient(correct_grad), - incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) { - if (nVar > MAXNVAR) { - SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION); - } - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar] = &JacobianBuffer[iVar * nVar]; - Jacobian_j[iVar] = &JacobianBuffer[iVar * nVar + MAXNVAR * MAXNVAR]; - } - - /*--- Initialize the JacobianBuffer to zero. ---*/ - for (unsigned short iVar = 0; iVar < 2*MAXNVAR*MAXNVAR; iVar++) { - JacobianBuffer[iVar] = 0.0; - } - } - - /*! - * \brief Compute the viscous residual using an average of gradients without correction. - * \param[in] config - Definition of the particular problem. - * \return A lightweight const-view (read-only) of the residual/flux and Jacobians. - */ - ResidualType<> ComputeResidual(const CConfig* config) final { - AD::StartPreacc(); - AD::SetPreaccIn(Coord_i, nDim); - AD::SetPreaccIn(Coord_j, nDim); - AD::SetPreaccIn(Normal, nDim); - AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); - AD::SetPreaccIn(ScalarVar_Grad_j, nVar, nDim); - if (correct_gradient) { - AD::SetPreaccIn(ScalarVar_i, nVar); - AD::SetPreaccIn(ScalarVar_j, nVar); - } - if (!std::is_same::value) { - AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]); - AD::SetPreaccIn(V_j[idx.Density()], V_j[idx.LaminarViscosity()], V_j[idx.EddyViscosity()]); - - Density_i = V_i[idx.Density()]; - Density_j = V_j[idx.Density()]; - Laminar_Viscosity_i = V_i[idx.LaminarViscosity()]; - Laminar_Viscosity_j = V_j[idx.LaminarViscosity()]; - Eddy_Viscosity_i = V_i[idx.EddyViscosity()]; - Eddy_Viscosity_j = V_j[idx.EddyViscosity()]; - } - - ExtraADPreaccIn(); - - su2double ProjGradScalarVarNoCorr[MAXNVAR]; - proj_vector_ij = ComputeProjectedGradient(nDim, nVar, Normal, Coord_i, Coord_j, ScalarVar_Grad_i, ScalarVar_Grad_j, - correct_gradient, ScalarVar_i, ScalarVar_j, ProjGradScalarVarNoCorr, - Proj_Mean_GradScalarVar); - FinishResidualCalc(config); - - AD::SetPreaccOut(Flux, nVar); - AD::EndPreacc(); - - return ResidualType<>(Flux, Jacobian_i, Jacobian_j); - } -}; diff --git a/SU2_CFD/include/numerics/scalar/scalar_edge_flux.hpp b/SU2_CFD/include/numerics/scalar/scalar_edge_flux.hpp new file mode 100644 index 00000000000..5753560e866 --- /dev/null +++ b/SU2_CFD/include/numerics/scalar/scalar_edge_flux.hpp @@ -0,0 +1,499 @@ +/*! + * \file scalar_edge_flux.hpp + * \brief Model-agnostic convection and diffusion of a transported scalar, shared by every + * scalar solver (turbulence, transition, species, flamelet, heat). + * \author P. Gomes + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, 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 . + */ + +#pragma once + +#include "../../../../Common/include/CConfig.hpp" +#include "../../../../Common/include/containers/container_decorators.hpp" +#include "../util.hpp" +#include "../../variables/CFlowVariable.hpp" + +/*! + * \brief Locates the data of one endpoint of an edge. + * \note The j side of a boundary flux indexes per-marker ghost containers, which have the + * same types as the solver's own, so the kernels read both through one code path. + */ +template +struct EdgeSide { + const VariableType& scalarNodes; /*!< \brief Scalar solver variables. */ + const CFlowVariable* flowNodes; /*!< \brief Flow variables, null for solid heat transfer. */ + CMatrixView coord; /*!< \brief Point coordinates. */ + CMatrixView gridVel; /*!< \brief Empty when the grid is static. */ +}; + +/*! + * \brief Loop invariant flags for a scalar edge flux, built once outside the edge loop so the + * compiler can unswitch the branches they guard. + * \note Built through the named constructors below: the flags are too many and too alike to be + * given positionally, where one transposed pair would change the discretization silently. + */ +struct ScalarFluxOptions { + bool dynamicGrid = false; + bool boundedScalar = false; + bool correctGradient = false; + bool accurateJacobians = false; + bool implicit = false; /*!< \brief Whether the Jacobians are assembled, and so computed. */ + bool convective = true; /*!< \brief Whether the convective scheme contributes. */ + bool viscous = false; /*!< \brief Whether the diffusion term contributes. */ + bool oneSided = false; /*!< \brief Whether only the row of i is assembled. */ + bool muscl = false; /*!< \brief Whether the convective scheme reconstructs. */ + + /*! + * \brief Options of the interior edge loop: both terms, both rows, and reconstruction and + * gradient correction as the configuration asks for them. + */ + static ScalarFluxOptions Interior(const CConfig& config, bool bounded, bool accurateJacobians = false) { + auto opt = common(config); + opt.boundedScalar = bounded; + opt.accurateJacobians = accurateJacobians; + opt.correctGradient = true; + opt.viscous = true; + opt.muscl = config.GetMUSCL(); + return opt; + } + + /*! + * \brief Options of a boundary that imposes a convective flux alone, which is most of them: + * the diffusive term at an inlet or an outlet causes serious convergence problems. + */ + static ScalarFluxOptions BoundaryConvective(const CConfig& config, bool bounded) { + auto opt = common(config); + opt.boundedScalar = bounded; + opt.oneSided = true; + return opt; + } + + /*! + * \brief Options of a boundary that imposes both terms, which is the turbomachinery sites. + * \note They impose no mass flux, so the bounded scheme contributes nothing here whatever the + * configuration says. + */ + static ScalarFluxOptions BoundaryFull(const CConfig& config) { + auto opt = common(config); + opt.correctGradient = true; + opt.viscous = true; + opt.oneSided = true; + return opt; + } + + /*! + * \brief Options of the diffusive pass of a fluid interface, which follows a convective pass + * over the donor vertices. + * \param[in] correctGrad - Whether the projected gradient is corrected for skewness, which the + * models do not agree on at this boundary. + */ + static ScalarFluxOptions BoundaryDiffusive(const CConfig& config, bool correctGrad) { + auto opt = common(config); + opt.correctGradient = correctGrad; + opt.convective = false; + opt.viscous = true; + opt.oneSided = true; + return opt; + } + + private: + static ScalarFluxOptions common(const CConfig& config) { + ScalarFluxOptions opt; + opt.dynamicGrid = config.GetDynamic_Grid(); + opt.implicit = config.GetKind_TimeIntScheme() == EULER_IMPLICIT; + return opt; + } +}; + +/*! + * \brief Thin wrapper giving a Vector the all(iVar) accessor the MUSCL reconstruction helpers + * expect (see CCompressiblePrimitives in numerics_simd/flow/variables.hpp). + */ +template +struct CScalarValues { + static constexpr size_t nVar = Size; /*!< \brief Only used as VarType::nVar when reconstruct's own nVarGrad_ default + applies; every call site here passes an explicit nVarGrad_ instead. */ + Vector all; +}; + +/*! + * \brief Diffusion of a transported scalar, driven by model-supplied coefficients. + * \note The derived class returns the coefficients of both orientations of the edge, as one + * object, so that whatever the two share is computed once. A model whose coefficients + * are the same in both orientations returns the same value twice. A model whose matrix + * is diagonal declares DiagonalDiffusion and returns a vector instead of a matrix. + */ +template +class CAvgGradScalarBase { + protected: + using Int = typename CLaneTraits::Int; + + /*! + * \param[in] rho - Density of both endpoints, read once by ComputeFlux. + */ + template + FORCEINLINE void diffusionTerms(const FlowIndices& idx, const ScalarFluxOptions& opt, Int iPoint, + const EdgeSide& side_i, Int jPoint, + const EdgeSide& side_j, const CPair& rho, + const Vector& normal, const Vector& vector_ij, + EdgeResidual& res) const { + if (!opt.viscous) return; + + constexpr size_t Size = EdgeResidual::Size; + + const Double dist2_ij = fmax(squaredNorm(vector_ij), EPS); + const Double proj_vector_ij = dot(vector_ij, normal) / dist2_ij; + + /*--- Average gradient, corrected for skewness when asked. + * \note Gathered one variable at a time, bounded by res.nVar rather than Size: a static + * model with nVar 1 has Size 2 (the Matrix degeneracy floor), and a + * dynamic one has Size MaxScalarVar, so a single Size-wide read would run past the + * actual width of the gradient container in either case. ---*/ + Matrix avgGrad; + for (size_t iVar = 0; iVar < res.nVar; ++iVar) { + const auto grad_i = gatherVariables(iPoint, side_i.scalarNodes.GetGradient(), iVar); + const auto grad_j = gatherVariables(jPoint, side_j.scalarNodes.GetGradient(), iVar); + for (int iDim = 0; iDim < nDim; ++iDim) avgGrad(iVar, iDim) = 0.5 * (grad_i(iDim) + grad_j(iDim)); + } + + if (opt.correctGradient) { + for (size_t iVar = 0; iVar < res.nVar; ++iVar) { + const Double phi_i = gatherVariables(iPoint, side_i.scalarNodes.GetSolution(), iVar); + const Double phi_j = gatherVariables(jPoint, side_j.scalarNodes.GetSolution(), iVar); + const Double corr = (dot(avgGrad[iVar], vector_ij) - phi_j + phi_i) / dist2_ij; + for (int iDim = 0; iDim < nDim; ++iDim) avgGrad(iVar, iDim) -= corr * vector_ij(iDim); + } + } + + Vector projGrad; + for (size_t iVar = 0; iVar < res.nVar; ++iVar) projGrad(iVar) = dot(avgGrad[iVar], normal); + + /*--- The Jacobians of a conservative model are w.r.t. the conserved (density-weighted) + * variable, which divides the geometric projection by the density of the row being written. ---*/ + Double proj_on_w_i = proj_vector_ij, proj_on_w_j = proj_vector_ij; + if constexpr (Derived::Conservative) { + proj_on_w_i = proj_vector_ij / rho.i; + proj_on_w_j = proj_vector_ij / rho.j; + } + + const auto* self = static_cast(this); + const auto D = self->coefficients(idx, iPoint, side_i, jPoint, side_j, rho); + + for (size_t iVar = 0; iVar < res.nVar; ++iVar) { + if constexpr (Derived::DiagonalDiffusion) { + res.flux_i(iVar) -= D.i(iVar) * projGrad(iVar); + if (opt.implicit) { + res.jac_ii(iVar, iVar) += D.i(iVar) * proj_on_w_i; + if (!opt.oneSided) res.jac_ij(iVar, iVar) -= D.i(iVar) * proj_on_w_j; + } + if (!opt.oneSided) { + res.flux_j(iVar) += D.j(iVar) * projGrad(iVar); + if (opt.implicit) { + res.jac_ji(iVar, iVar) -= D.j(iVar) * proj_on_w_i; + res.jac_jj(iVar, iVar) += D.j(iVar) * proj_on_w_j; + } + } + } else { + for (size_t jVar = 0; jVar < res.nVar; ++jVar) { + res.flux_i(iVar) -= D.i(iVar, jVar) * projGrad(jVar); + if (opt.implicit) { + res.jac_ii(iVar, jVar) += D.i(iVar, jVar) * proj_on_w_i; + if (!opt.oneSided) res.jac_ij(iVar, jVar) -= D.i(iVar, jVar) * proj_on_w_j; + } + if (!opt.oneSided) { + res.flux_j(iVar) += D.j(iVar, jVar) * projGrad(jVar); + if (opt.implicit) { + res.jac_ji(iVar, jVar) -= D.j(iVar, jVar) * proj_on_w_i; + res.jac_jj(iVar, jVar) += D.j(iVar, jVar) * proj_on_w_j; + } + } + } + } + } + + if (opt.implicit && opt.accurateJacobians) { + /*--- Coefficients that depend on the transported variables contribute here, from whatever + * the model chose to carry in the object it returned from coefficients. ---*/ + self->coefficientJacobians(opt, D, projGrad, res); + } + + self->extraDiffusionTerms(idx, opt, iPoint, side_i, jPoint, side_j, rho, normal, vector_ij, res); + } + + /*! + * \brief Contribution of the derivatives of the coefficients themselves. + */ + template + FORCEINLINE void coefficientJacobians(Ts&&...) const {} + + /*! + * \brief Diffusion of a model that transports more than one gradient, of states it + * synthesises from its own containers. + */ + template + FORCEINLINE void extraDiffusionTerms(Ts&&...) const {} +}; + +/*! + * \brief Convective flux shared by every model whose transport equation has the shape + * flux(iVar) = a0 * w_i * phi_i(iVar) + a1 * w_j * phi_j(iVar), which is every + * model except SA and stochastic backscatter (see CScalarFlux_SA). + * \note The weight w is 1 for a non-conservative model, the density for a conservative one. + */ +template +class CUpwScalarFlux : public CAvgGradScalarBase { + protected: + using Int = typename CLaneTraits::Int; + + explicit CUpwScalarFlux(const CConfig&) {} + + /*! + * \brief Upwind convection of the transported variable, weighted by the density for a + * conservative model. + * \param[in] phi - Transported variable of both endpoints, reconstructed if opt.muscl is set; + * read from here rather than side_i/side_j.scalarNodes directly so a model needs + * no reconstruction logic of its own. + * \param[in] rho - Density of both endpoints, reconstructed alongside the velocity when the + * convective scheme reconstructs, so it weights the flux as the velocity does. + * \note The flux is written in terms of the transported variable but the Jacobians are w.r.t. + * the conserved one, which for a conservative model is the density-weighted variable; + * the density therefore multiplies the flux and not the Jacobian. + */ + template + FORCEINLINE void finalizeFlux(const FlowIndices&, const ScalarFluxOptions& opt, Int, const EdgeSide&, + Int, const EdgeSide&, const Double& a0, const Double& a1, + const CPair& rho, const CPair>& phi, + EdgeResidual& res) const { + Double w0 = a0, w1 = a1; + if constexpr (Derived::Conservative) { + w0 *= rho.i; + w1 *= rho.j; + } + + for (size_t iVar = 0; iVar < res.nVar; ++iVar) { + const Double flux = w0 * phi.i.all(iVar) + w1 * phi.j.all(iVar); + + res.flux_i(iVar) += flux; + if (!opt.oneSided) res.flux_j(iVar) -= flux; + + if (opt.implicit) { + res.jac_ii(iVar, iVar) += a0; + if (!opt.oneSided) { + res.jac_ij(iVar, iVar) += a1; + res.jac_ji(iVar, iVar) -= a0; + res.jac_jj(iVar, iVar) -= a1; + } + } + } + } +}; + +/*! + * \brief Upwind convection and diffusion of a transported scalar, accumulated into one + * residual, each term contributing or not according to the options. + */ +template +class CUpwScalarBase : public CUpwScalarFlux { + public: + using Double = Double_; + using Int = typename CLaneTraits::Int; + static constexpr int nDim = nDim_; + static constexpr size_t nVar = nVar_; + + /*! + * \brief Backing size, for the model to size the containers its coefficients() returns; + * nVar itself is Dynamic for a runtime model and never usable as a container size. + */ + static constexpr size_t Size = EdgeResidual::Size; + + /*! + * \brief Whether the model's diffusion coefficients read the density, beyond the reading that + * Conservative already implies. A model that declares neither never gathers it. + */ + static constexpr bool DiffusionReadsDensity = false; + + protected: + using Base = CUpwScalarFlux; + + const FlowIndices idx; + const size_t nEqn; /*!< \brief Equations of the model, which a dynamic one gives to its base. */ + + /*! + * \brief MUSCL reconstruction parameters, read from CConfig once per construction (i.e. once + * per nonlinear iteration, see CScalarSolver::EdgeFluxResidual) instead of per edge; + * this is also where the scalar limiter's freezing (GetLimiterIter) is resolved, by + * collapsing its type to NONE once frozen. The flow limiter is not frozen this way: once + * the flow solver stops recomputing it, it keeps applying the last values it has. + */ + const su2double kappa, umusclRamp, kappaFlow; + const LIMITER limiterType, limiterTypeFlow; + const bool musclFlow; + + public: + /*! + * \brief Constructor, inherited by the model with `using Base::Base`. + * \note Public, not protected: a using-declaration that inherits a constructor keeps the + * base's own access, so the solver that builds the concrete model needs this public + * to build it at all. + */ + explicit CUpwScalarBase(const CConfig& config, size_t nEqn_ = nVar_) + : Base(config), + idx(nDim, config.GetnSpecies()), + nEqn(nEqn_), + kappa(config.GetMUSCL_Kappa()), + umusclRamp(config.GetMUSCLRampValue()), + kappaFlow(config.GetMUSCL_Kappa_Flow()), + limiterType(config.GetInnerIter() <= config.GetLimiterIter() ? config.GetKind_SlopeLimit() : LIMITER::NONE), + limiterTypeFlow(config.GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE ? config.GetKind_SlopeLimit_Flow() + : LIMITER::NONE), + musclFlow(config.GetMUSCL_Flow() && config.GetKind_ConvNumScheme_Flow() == SPACE_UPWIND) { + if (nEqn > Size) { + SU2_MPI::Error("Static arrays are too small for the requested equation count.", CURRENT_FUNCTION); + } + } + + template + FORCEINLINE EdgeResidual ComputeFlux(const ScalarFluxOptions& opt, Int iPoint, + const EdgeSide& side_i, Int jPoint, + const EdgeSide& side_j, + const Vector& normal, const Double& massFlux) const { + /*--- Inputs are registered as they are read, by each of the two terms. ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(normal, nDim); + + EdgeResidual res(nEqn); + + /*--- Read once by the reconstruction and by the diffusion. ---*/ + Vector vector_ij; + if (opt.muscl || opt.viscous) { + vector_ij = distanceVector(iPoint, side_i.coord, jPoint, side_j.coord); + } + + /*--- Density of both endpoints, gathered once: the conservative weighting of the convective + * term, the bounded scheme's division of the mass flux, and some models' diffusion + * coefficients all want it, and in reverse mode every gather is a preaccumulation input. ---*/ + CPair rho{Double(1.0), Double(1.0)}; + if (Derived::Conservative || Derived::DiffusionReadsDensity || opt.boundedScalar) { + rho.i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.Density()); + rho.j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.Density()); + } + + if (opt.convective) { + /*--- Upwinding weights of the face normal mass or volume flux, and the density that weights + * a conservative flux, which follows the velocity in being reconstructed or not. ---*/ + Double a0, a1; + CPair rhoConv = rho; + + if (opt.boundedScalar) { + AD::SetPreaccIn(massFlux); + a0 = fmax(0.0, massFlux) / rho.i; + a1 = fmin(0.0, massFlux) / rho.j; + } else { + CPair> u; + u.i.all = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.Velocity()); + u.j.all = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.Velocity()); + + if (opt.muscl && musclFlow) { + reconstruct(iPoint, jPoint, vector_ij, side_i.flowNodes->GetGradient_Reconstruction(), + side_i.flowNodes->GetLimiter_Primitive(), limiterTypeFlow, idx.Velocity(), u, kappaFlow, + umusclRamp); + + if constexpr (Derived::Conservative) { + /*--- Density is not adjacent to the velocity in the primitive row, so it is a second + * reconstruction of one variable rather than a wider slice of the first. ---*/ + CPair> r; + r.i.all(0) = rho.i; + r.j.all(0) = rho.j; + reconstruct<1>(iPoint, jPoint, vector_ij, side_i.flowNodes->GetGradient_Reconstruction(), + side_i.flowNodes->GetLimiter_Primitive(), limiterTypeFlow, idx.Density(), r, kappaFlow, + umusclRamp); + rhoConv.i = r.i.all(0); + rhoConv.j = r.j.all(0); + } + } + + /*--- Face normal velocity of the mean of the two points, relative to the grid. ---*/ + Vector vel_ij; + for (int iDim = 0; iDim < nDim; ++iDim) vel_ij(iDim) = 0.5 * (u.i.all(iDim) + u.j.all(iDim)); + + if (opt.dynamicGrid) { + const auto ug_i = gatherVariables(iPoint, side_i.gridVel); + /*--- A boundary's ghost point has no grid velocity of its own: it is spatially + * coincident with i, so it moves with it. ---*/ + const auto ug_j = opt.oneSided ? ug_i : gatherVariables(jPoint, side_j.gridVel); + for (int iDim = 0; iDim < nDim; ++iDim) vel_ij(iDim) -= 0.5 * (ug_i(iDim) + ug_j(iDim)); + } + + const Double q_ij = dot(vel_ij, normal); + a0 = fmax(0.0, q_ij); + a1 = fmin(0.0, q_ij); + } + + /*--- Transported variable of both endpoints, reconstructed if opt.muscl is set. + * Gathered one variable at a time (like the diffusion gradients above) so a static model + * with nVar 1 never reads past the single column its solution container actually has. ---*/ + CPair> phi; + for (size_t iVar = 0; iVar < res.nVar; ++iVar) { + phi.i.all(iVar) = gatherVariables(iPoint, side_i.scalarNodes.GetSolution(), iVar); + phi.j.all(iVar) = gatherVariables(jPoint, side_j.scalarNodes.GetSolution(), iVar); + } + + if (opt.muscl) { + if constexpr (nVar != Dynamic) { + reconstruct(iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(), + side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp); + } else { + /*--- A dynamic model's equation count is only known at runtime, so the reconstructed + * width is passed as an argument instead of a template parameter. ---*/ + reconstruct(iPoint, jPoint, vector_ij, side_i.scalarNodes.GetGradient_Reconstruction(), + side_i.scalarNodes.GetLimiter(), limiterType, 0, phi, kappa, umusclRamp, res.nVar); + } + } + + static_cast(this)->finalizeFlux(idx, opt, iPoint, side_i, jPoint, side_j, a0, a1, rhoConv, phi, + res); + } + + Base::diffusionTerms(idx, opt, iPoint, side_i, jPoint, side_j, rho, normal, vector_ij, res); + + setPreaccOut(res.flux_i, res.nVar); + if (!opt.oneSided) setPreaccOut(res.flux_j, res.nVar); + AD::EndPreacc(); + + return res; + } + + /*! + * \brief Compute the flux of an edge and write it to the linear system. + */ + template + FORCEINLINE void ComputeFlux(const ScalarFluxOptions& opt, Int iEdge, Int iPoint, + const EdgeSide& side_i, Int jPoint, const EdgeSide& side_j, + const Vector& normal, const Double& massFlux, UpdateType updateType, + Double updateMask, CSysVector& vector, CSysVector& vectorDiff, + SparseMatrixType& matrix) const { + const auto res = ComputeFlux(opt, iPoint, side_i, jPoint, side_j, normal, massFlux); + + updateLinearSystem(iEdge, iPoint, jPoint, opt.implicit, updateType, updateMask, res, vector, vectorDiff, matrix); + } +}; diff --git a/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp b/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp new file mode 100644 index 00000000000..dc634a235ca --- /dev/null +++ b/SU2_CFD/include/numerics/species/flamelet_edge_flux.hpp @@ -0,0 +1,172 @@ +/*! + * \file flamelet_edge_flux.hpp + * \brief Flamelet transport as a third-layer scalar flux, see numerics/scalar/scalar_edge_flux.hpp. + * \author P. Gomes + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, 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 . + */ + +#pragma once + +#include "species_edge_flux.hpp" + +/*! + * \class CScalarFlux_Flamelet + * \ingroup ViscDiscr + * \brief Convection and diffusion of the flamelet controlling variables and passive species, + * which is species transport plus two preferential diffusion terms. + * \note The preferential diffusion terms read the beta scalars and their gradients from the + * auxiliary variables of the solver's own containers, which the per-marker ghost containers + * of a boundary do not carry. They are an interior edge term: boundaries instantiate + * CScalarFlux_Species, as they did before this model existed. + */ +template +class CScalarFlux_Flamelet final + : public CScalarFluxSpeciesBase, FlowIndices, nDim, + nVar> { + public: + using Base = CScalarFluxSpeciesBase; + using Int = typename Base::Int; + + explicit CScalarFlux_Flamelet(const CConfig& config) + : Base(config), + preferentialDiffusion(config.GetFlameletParsedOptions().preferential_diffusion), + nControlVars(config.GetFlameletParsedOptions().n_control_vars) {} + + /*! + * \brief Preferential diffusion, two terms with the shape of the ordinary diffusion but of + * states the model synthesises: div(D grad(beta - phi)) for each controlling variable, + * and a thermal term div(beta_T D grad(T)) on the enthalpy equation. + * \note The thermal term has no implicit part, matching the treatment of the heat flux it + * models; the first term has the same thin shear layer Jacobian as the ordinary diffusion, + * because it is the same operator applied to a shifted state. + */ + template + FORCEINLINE void extraDiffusionTerms(const FlowIndices& idx, const ScalarFluxOptions& opt, Int iPoint, + const EdgeSide& side_i, Int jPoint, + const EdgeSide& side_j, const CPair& rho, + const Vector& normal, const Vector& vector_ij, + EdgeResidual& res) const { + if (!preferentialDiffusion) return; + + const Double dist2_ij = fmax(squaredNorm(vector_ij), EPS); + const Double proj_vector_ij = dot(vector_ij, normal) / dist2_ij; + const Double proj_on_rho_i = proj_vector_ij / rho.i; + const Double proj_on_rho_j = proj_vector_ij / rho.j; + + const Double diffTurb = Base::turbulentDiffusivity(idx, iPoint, side_i, jPoint, side_j); + + /*--- The gradient of a controlling variable is subtracted from that of its beta scalar, so + * that what is added here is the difference from the ordinary diffusion already applied. ---*/ + for (auto iScalar = 0u; iScalar < nControlVars; ++iScalar) { + const auto iBeta = betaIndex(iScalar); + + const Double phi_i = gatherVariables(iPoint, side_i.scalarNodes.GetAuxVar(), iBeta) - + gatherVariables(iPoint, side_i.scalarNodes.GetSolution(), iScalar); + const Double phi_j = gatherVariables(jPoint, side_j.scalarNodes.GetAuxVar(), iBeta) - + gatherVariables(jPoint, side_j.scalarNodes.GetSolution(), iScalar); + + auto grad_i = gatherVariables(iPoint, side_i.scalarNodes.GetAuxVarGradient(), iBeta); + auto grad_j = gatherVariables(jPoint, side_j.scalarNodes.GetAuxVarGradient(), iBeta); + const auto gradPhi_i = gatherVariables(iPoint, side_i.scalarNodes.GetGradient(), iScalar); + const auto gradPhi_j = gatherVariables(jPoint, side_j.scalarNodes.GetGradient(), iScalar); + for (int iDim = 0; iDim < nDim; ++iDim) { + grad_i(iDim) -= gradPhi_i(iDim); + grad_j(iDim) -= gradPhi_j(iDim); + } + + const Double D_i = gatherVariables(iPoint, side_i.scalarNodes.GetDiffusivity(), iScalar); + const Double D_j = gatherVariables(jPoint, side_j.scalarNodes.GetDiffusivity(), iScalar); + const Double D = 0.5 * (rho.i * D_i + rho.j * D_j) + diffTurb; + + const Double projGrad = projectedGradient(opt, grad_i, grad_j, phi_i, phi_j, normal, vector_ij, dist2_ij); + + res.flux_i(iScalar) -= D * projGrad; + if (!opt.oneSided) res.flux_j(iScalar) += D * projGrad; + + if (opt.implicit) { + res.jac_ii(iScalar, iScalar) += D * proj_on_rho_i; + if (!opt.oneSided) { + res.jac_ij(iScalar, iScalar) -= D * proj_on_rho_j; + res.jac_ji(iScalar, iScalar) -= D * proj_on_rho_i; + res.jac_jj(iScalar, iScalar) += D * proj_on_rho_j; + } + } + } + + /*--- Thermal term, on the enthalpy equation alone, driven by the temperature gradient. ---*/ + if (nControlVars <= I_ENTH) return; + + const Double T_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.Temperature()); + const Double T_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.Temperature()); + + const auto gradT_i = gatherVariables(iPoint, side_i.flowNodes->GetGradient_Primitive(), idx.Temperature()); + const auto gradT_j = gatherVariables(jPoint, side_j.flowNodes->GetGradient_Primitive(), idx.Temperature()); + + const Double Dth_i = gatherVariables(iPoint, side_i.scalarNodes.GetAuxVar(), I_BETA_ENTH_THERMAL) * + gatherVariables(iPoint, side_i.scalarNodes.GetDiffusivity(), I_ENTH); + const Double Dth_j = gatherVariables(jPoint, side_j.scalarNodes.GetAuxVar(), I_BETA_ENTH_THERMAL) * + gatherVariables(jPoint, side_j.scalarNodes.GetDiffusivity(), I_ENTH); + const Double Dth = 0.5 * (rho.i * Dth_i + rho.j * Dth_j) + diffTurb; + + const Double projGradT = projectedGradient(opt, gradT_i, gradT_j, T_i, T_j, normal, vector_ij, dist2_ij); + + res.flux_i(I_ENTH) -= Dth * projGradT; + if (!opt.oneSided) res.flux_j(I_ENTH) += Dth * projGradT; + } + + private: + const bool preferentialDiffusion; + const unsigned short nControlVars; + + /*! + * \brief Auxiliary variable holding the beta scalar of a controlling variable. + */ + static FORCEINLINE unsigned short betaIndex(unsigned short iScalar) { + switch (iScalar) { + case I_PROGVAR: + return I_BETA_PROGVAR; + case I_ENTH: + return I_BETA_ENTH; + default: + return I_BETA_MIXFRAC; + } + } + + /*! + * \brief Average gradient of one synthesised state projected on the normal, corrected for + * skewness when asked, which is what the ordinary diffusion does for a transported one. + */ + FORCEINLINE Double projectedGradient(const ScalarFluxOptions& opt, const Vector& grad_i, + const Vector& grad_j, const Double& phi_i, const Double& phi_j, + const Vector& normal, const Vector& vector_ij, + const Double& dist2_ij) const { + Vector avgGrad; + for (int iDim = 0; iDim < nDim; ++iDim) avgGrad(iDim) = 0.5 * (grad_i(iDim) + grad_j(iDim)); + + if (opt.correctGradient) { + const Double corr = (dot(avgGrad, vector_ij) - phi_j + phi_i) / dist2_ij; + for (int iDim = 0; iDim < nDim; ++iDim) avgGrad(iDim) -= corr * vector_ij(iDim); + } + return dot(avgGrad, normal); + } +}; diff --git a/SU2_CFD/include/numerics/species/species_convection.hpp b/SU2_CFD/include/numerics/species/species_convection.hpp deleted file mode 100644 index 114501cb544..00000000000 --- a/SU2_CFD/include/numerics/species/species_convection.hpp +++ /dev/null @@ -1,85 +0,0 @@ -/*! - * \file species_convection.hpp - * \brief Declarations of numerics classes for discretization of - * convective fluxes in species problems. - * \author T. Kattmann - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, 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 . - */ - -#pragma once - -#include "../scalar/scalar_convection.hpp" - -/*! - * \class CUpwSca_Species - * \brief Class for doing a scalar upwind solver for the species transport equations. - * \ingroup ConvDiscr - */ -template -class CUpwSca_Species final : public CUpwScalar { - private: - using Base = CUpwScalar; - using Base::nVar; - using Base::nDim; - using Base::V_i; - using Base::V_j; - using Base::a0; - using Base::a1; - using Base::Flux; - using Base::Jacobian_i; - using Base::Jacobian_j; - using Base::ScalarVar_i; - using Base::ScalarVar_j; - using Base::idx; - using Base::bounded_scalar; - - /*! - * \brief Adds any extra variables to AD - */ - void ExtraADPreaccIn() override {} - - /*! - * \brief Species transport specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - for (auto iVar = 0u; iVar < nVar; iVar++) { - Flux[iVar] = a0 * V_i[idx.Density()] * ScalarVar_i[iVar] + a1 * V_j[idx.Density()] * ScalarVar_j[iVar]; - - /*--- Jacobians are taken wrt rho*Y not Y alone in the species solver. ---*/ - /*--- Off-diagonal entries are zero. ---*/ - Jacobian_i[iVar][iVar] = a0; - Jacobian_j[iVar][iVar] = a1; - } // iVar - } - - public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CUpwSca_Species(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) - : CUpwScalar(val_nDim, val_nVar, config) { bounded_scalar = config->GetBounded_Species(); } -}; diff --git a/SU2_CFD/include/numerics/species/species_diffusion.hpp b/SU2_CFD/include/numerics/species/species_diffusion.hpp deleted file mode 100644 index 424c4a47b8c..00000000000 --- a/SU2_CFD/include/numerics/species/species_diffusion.hpp +++ /dev/null @@ -1,110 +0,0 @@ -/*! - * \file species_diffusion.hpp - * \brief Declarations of numerics classes for discretization of - * viscous fluxes in species problems. - * \author T. Kattmann - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, 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 . - */ - -#pragma once - -#include "../scalar/scalar_diffusion.hpp" - -/*! - * \class CAvgGrad_Species - * \brief Class for computing viscous term using average of gradients (species transport model). - * \ingroup ViscDiscr - */ -template -class CAvgGrad_Species final : public CAvgGrad_Scalar { - private: - using Base = CAvgGrad_Scalar; - using Base::nVar; - using Base::Eddy_Viscosity_i; - using Base::Eddy_Viscosity_j; - using Base::Diffusion_Coeff_i; - using Base::Diffusion_Coeff_j; - using Base::Density_i; - using Base::Density_j; - using Base::ScalarVar_i; - using Base::ScalarVar_j; - using Base::Proj_Mean_GradScalarVar; - using Base::proj_vector_ij; - using Base::Flux; - using Base::Jacobian_i; - using Base::Jacobian_j; - - const bool turbulence; - - /*! - * \brief Adds any extra variables to AD - */ - void ExtraADPreaccIn(void) override { - AD::SetPreaccIn(Diffusion_Coeff_i, nVar); - AD::SetPreaccIn(Diffusion_Coeff_j, nVar); - } - - /*! - * \brief Species transport specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - for (auto iVar = 0u; iVar < nVar; iVar++) { - - const su2double Diffusivity_Lam = 0.5 * (Density_i * Diffusion_Coeff_i[iVar] + Density_j * Diffusion_Coeff_j[iVar]); - - su2double Diffusivity_Turb = 0.0; - - if (turbulence) { - const su2double Sc_t = config->GetSchmidt_Number_Turbulent(); - Diffusivity_Turb = 0.5 * (Eddy_Viscosity_i / Sc_t + Eddy_Viscosity_j / Sc_t); - } - - const su2double Diffusivity = Diffusivity_Lam + Diffusivity_Turb; - - Flux[iVar] = Diffusivity * Proj_Mean_GradScalarVar[iVar]; - - /*--- Use TSL approx. to compute derivatives of the gradients. ---*/ - - /*--- Off-diagonal entries are all zero. ---*/ - const su2double proj_on_rhoi = proj_vector_ij / Density_i; - Jacobian_i[iVar][iVar] = -Diffusivity * proj_on_rhoi; - - const su2double proj_on_rhoj = proj_vector_ij / Density_j; - Jacobian_j[iVar][iVar] = Diffusivity * proj_on_rhoj; - - } // iVar - } - - public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] correct_grad - Whether to correct gradient for skewness. - * \param[in] config - Definition of the particular problem. - */ - CAvgGrad_Species(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config) - : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config), - turbulence(config->GetKind_Turb_Model() != TURB_MODEL::NONE) {} -}; diff --git a/SU2_CFD/include/numerics/species/species_edge_flux.hpp b/SU2_CFD/include/numerics/species/species_edge_flux.hpp new file mode 100644 index 00000000000..e744d8c9655 --- /dev/null +++ b/SU2_CFD/include/numerics/species/species_edge_flux.hpp @@ -0,0 +1,113 @@ +/*! + * \file species_edge_flux.hpp + * \brief Species transport as a third-layer scalar flux, see numerics/scalar/scalar_edge_flux.hpp. + * \author P. Gomes + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, 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 . + */ + +#pragma once + +#include "../scalar/scalar_edge_flux.hpp" + +/*! + * \class CScalarFluxSpeciesBase + * \ingroup ViscDiscr + * \brief Convection and diffusion of a mass fraction, conservative with a diagonal, + * i/j-symmetric diffusion matrix. Unlike SA/SST/LM, the equation count is only known at + * runtime (one per transported species), so this is the framework's first Dynamic-nVar + * model: nEqn is passed to the base explicitly, and coefficients() loops to it rather than + * to a compile-time nVar. + * \note This carries no finalizeFlux of its own: the inherited CUpwScalarFlux one is exactly + * flux(iVar) = a0*rho_i*Y_i(iVar) + a1*rho_j*Y_j(iVar), Conservative weighting by density, + * which is the model's whole convective term. + * \note It takes the most derived class as a parameter so that the flamelet model, which adds a + * preferential diffusion term to the same coefficients, is a sibling rather than a copy. + */ +template +class CScalarFluxSpeciesBase : public CUpwScalarBase { + public: + static constexpr bool Conservative = true; + static constexpr bool DiagonalDiffusion = true; + + using Base = CUpwScalarBase; + using Int = typename Base::Int; + + explicit CScalarFluxSpeciesBase(const CConfig& config) + : Base(config, config.GetnSpecies()), + turbulence(config.GetKind_Turb_Model() != TURB_MODEL::NONE), + Sc_t(config.GetSchmidt_Number_Turbulent()) {} + + /*! + * \brief Diffusion coefficients, an i/j average of (rho * mass diffusivity) per species, plus a + * turbulent (mu_t/Sc_t) contribution shared by every species, when a turbulence model is + * active; identical for both edge sides. + */ + template + FORCEINLINE CPair> coefficients(const FlowIndices& idx, Int iPoint, + const EdgeSide& side_i, Int jPoint, + const EdgeSide& side_j, + const CPair& rho) const { + const Double diffTurb = turbulentDiffusivity(idx, iPoint, side_i, jPoint, side_j); + + Vector D; + for (size_t iVar = 0; iVar < this->nEqn; ++iVar) { + const Double D_lam_i = gatherVariables(iPoint, side_i.scalarNodes.GetDiffusivity(), iVar); + const Double D_lam_j = gatherVariables(jPoint, side_j.scalarNodes.GetDiffusivity(), iVar); + D(iVar) = 0.5 * (rho.i * D_lam_i + rho.j * D_lam_j) + diffTurb; + } + return {D, D}; + } + + protected: + /*! + * \brief Turbulent contribution to the diffusivity, shared by every species and, in the + * flamelet model, by the preferential diffusion terms. + */ + template + FORCEINLINE Double turbulentDiffusivity(const FlowIndices& idx, Int iPoint, const EdgeSide& side_i, + Int jPoint, const EdgeSide& side_j) const { + if (!turbulence) return Double(0.0); + + const Double muT_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.EddyViscosity()); + const Double muT_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.EddyViscosity()); + return 0.5 * (muT_i / Sc_t + muT_j / Sc_t); + } + + private: + const bool turbulence; + const su2double Sc_t; +}; + +/*! + * \class CScalarFlux_Species + * \ingroup ViscDiscr + * \brief Convection and diffusion of the species transport model. + */ +template +class CScalarFlux_Species final + : public CScalarFluxSpeciesBase, FlowIndices, nDim, + nVar> { + public: + using Base = CScalarFluxSpeciesBase; + using Base::Base; +}; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp deleted file mode 100644 index 2338c2ca837..00000000000 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_convection.hpp +++ /dev/null @@ -1,40 +0,0 @@ -/*! - * \file trans_convection.hpp - * \brief Delarations of numerics classes for discretization of - * convective fluxes in transition problems. - * \author S. Kang - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, 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 . - */ - -#pragma once - -#include "../turb_convection.hpp" - -/*! - * \class CUpwSca_TransLM - * \brief Re-use the SST convective fluxes for the scalar upwind discretization of LM transition model equations. - * \ingroup ConvDiscr - */ -template -using CUpwSca_TransLM = CUpwSca_TurbSST; - diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp deleted file mode 100644 index 0d4aee5a947..00000000000 --- a/SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp +++ /dev/null @@ -1,105 +0,0 @@ -/*! - * \file trans_diffusion.hpp - * \brief Declarations of numerics classes for discretization of - * viscous fluxes in transition problems. - * \author S. Kang - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, 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 . - */ -#pragma once - - -#include "../../scalar/scalar_diffusion.hpp" - -/*! - * \class CAvgGrad_TransLM - * \brief Class for computing viscous term using average of gradient with correction (LM transition model). - * \ingroup ViscDiscr - * \author S. Kang. - */ -template -class CAvgGrad_TransLM final : public CAvgGrad_Scalar { -private: - using Base = CAvgGrad_Scalar; - using Base::Laminar_Viscosity_i; - using Base::Laminar_Viscosity_j; - using Base::Eddy_Viscosity_i; - using Base::Eddy_Viscosity_j; - using Base::Density_i; - using Base::Density_j; - using Base::ScalarVar_i; - using Base::ScalarVar_j; - using Base::Proj_Mean_GradScalarVar; - using Base::proj_vector_ij; - using Base::Flux; - using Base::Jacobian_i; - using Base::Jacobian_j; - - /*! - * \brief Adds any extra variables to AD - */ - void ExtraADPreaccIn() override {} - - /*! - * \brief LM transition model specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; - - /*--- Compute mean effective dynamic viscosity ---*/ - const su2double diff_i_gamma = Laminar_Viscosity_i + Eddy_Viscosity_i; - const su2double diff_j_gamma = Laminar_Viscosity_j + Eddy_Viscosity_j; - const su2double diff_i_ReThetaT = 2.0*(Laminar_Viscosity_i + Eddy_Viscosity_i); - const su2double diff_j_ReThetaT = 2.0*(Laminar_Viscosity_j + Eddy_Viscosity_j); - - const su2double diff_gamma = 0.5*(diff_i_gamma + diff_j_gamma); - const su2double diff_ReThetaT = 0.5*(diff_i_ReThetaT + diff_j_ReThetaT); - - Flux[0] = diff_gamma*Proj_Mean_GradScalarVar[0]; - Flux[1] = diff_ReThetaT*Proj_Mean_GradScalarVar[1]; - - /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/ - if (implicit) { - const su2double proj_on_rho_i = proj_vector_ij/Density_i; - Jacobian_i[0][0] = -diff_gamma*proj_on_rho_i; Jacobian_i[0][1] = 0.0; - Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_ReThetaT*proj_on_rho_i; - - const su2double proj_on_rho_j = proj_vector_ij/Density_j; - Jacobian_j[0][0] = diff_gamma*proj_on_rho_j; Jacobian_j[0][1] = 0.0; - Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_ReThetaT*proj_on_rho_j; - } - } - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] correct_grad - Whether to correct gradient for skewness. - * \param[in] config - Definition of the particular problem. - */ - CAvgGrad_TransLM(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config) - : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config){ - } - -}; diff --git a/SU2_CFD/include/numerics/turbulent/transition/trans_edge_flux.hpp b/SU2_CFD/include/numerics/turbulent/transition/trans_edge_flux.hpp new file mode 100644 index 00000000000..927fcb6f5b0 --- /dev/null +++ b/SU2_CFD/include/numerics/turbulent/transition/trans_edge_flux.hpp @@ -0,0 +1,78 @@ +/*! + * \file trans_edge_flux.hpp + * \brief Langtry-Menter transition model as a third-layer scalar flux, see numerics/scalar/scalar_edge_flux.hpp. + * \author P. Gomes + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, 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 . + */ + +#pragma once + +#include "../../scalar/scalar_edge_flux.hpp" + +/*! + * \class CScalarFlux_TransLM + * \ingroup ViscDiscr + * \brief Convection and diffusion of the Langtry-Menter transition model, conservative with a + * diagonal, i/j-symmetric diffusion matrix. The coefficients depend only on the flow's + * mu/mu_t, not on the transported gamma/Re_theta, so no coefficientJacobians override + * is needed. + * \note LM writes no finalizeFlux of its own: the inherited CUpwScalarFlux one is exactly + * flux(iVar) = a0*rho_i*phi_i(iVar) + a1*rho_j*phi_j(iVar), Conservative weighting by + * density, which is the model's whole convective term. + */ +template +class CScalarFlux_TransLM + : public CUpwScalarBase, FlowIndices, nDim, nVar> { + public: + static constexpr bool Conservative = true; + static constexpr bool DiagonalDiffusion = true; + + using Base = CUpwScalarBase; + using Int = typename Base::Int; + using Base::Base; + + /*! + * \brief Diffusion coefficients, an i/j average of (mu+mu_t) for intermittency and of + * 2*(mu+mu_t) for the momentum-thickness Reynolds number; identical for both edge sides. + */ + template + FORCEINLINE CPair> coefficients(const FlowIndices& idx, Int iPoint, + const EdgeSide& side_i, Int jPoint, + const EdgeSide& side_j, + const CPair&) const { + const Double mu_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.LaminarViscosity()); + const Double mu_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.LaminarViscosity()); + const Double muT_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.EddyViscosity()); + const Double muT_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.EddyViscosity()); + + const Double diff_i_gamma = mu_i + muT_i; + const Double diff_j_gamma = mu_j + muT_j; + const Double diff_i_ReThetaT = 2.0 * (mu_i + muT_i); + const Double diff_j_ReThetaT = 2.0 * (mu_j + muT_j); + + Vector D; + D(0) = 0.5 * (diff_i_gamma + diff_j_gamma); + D(1) = 0.5 * (diff_i_ReThetaT + diff_j_ReThetaT); + return {D, D}; + } +}; diff --git a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp deleted file mode 100644 index 6c1641db87d..00000000000 --- a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp +++ /dev/null @@ -1,141 +0,0 @@ -/*! - * \file turb_convection.hpp - * \brief Declarations of numerics classes for discretization of - * convective fluxes in turbulence problems. - * \author F. Palacios, T. Economon - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, 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 . - */ - -#pragma once - -#include "../scalar/scalar_convection.hpp" - -/*! - * \class CUpwSca_TurbSA - * \brief Class for doing a scalar upwind solver for the Spalar-Allmaras turbulence model equations. - * \ingroup ConvDiscr - * \author A. Bueno. - */ -template -class CUpwSca_TurbSA final : public CUpwScalar { -private: - using Base = CUpwScalar; - using Base::a0; - using Base::a1; - using Base::Flux; - using Base::Jacobian_i; - using Base::Jacobian_j; - using Base::ScalarVar_i; - using Base::ScalarVar_j; - using Base::bounded_scalar; - using Base::V_i; - using Base::V_j; - using Base::idx; - using Base::nVar; - - /*! - * \brief Adds any extra variables to AD. - */ - void ExtraADPreaccIn() override {} - - /*! - * \brief SA specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - if (config->GetSBSParam().StochasticBackscatter && config->GetSBSParam().SBS_Ctau > 0.0) { - for (unsigned short iVar = 1; iVar < nVar; iVar++) { - Flux[iVar] = (a0 + a1) * 0.5 * (ScalarVar_i[iVar] + ScalarVar_j[iVar]); - Jacobian_i[iVar][iVar] = 0.5 * (a0+a1); - Jacobian_j[iVar][iVar] = 0.5 * (a0+a1); - } - } - Flux[0] = a0*ScalarVar_i[0] + a1*ScalarVar_j[0]; - Jacobian_i[0][0] = a0; - Jacobian_j[0][0] = a1; - } - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CUpwSca_TurbSA(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) - : CUpwScalar(val_nDim, val_nVar, config) { bounded_scalar = config->GetBounded_Turb(); } -}; - -/*! - * \class CUpwSca_TurbSST - * \brief Class for doing a scalar upwind solver for the Menter SST turbulence model equations. - * \ingroup ConvDiscr - * \author A. Campos. - */ -template -class CUpwSca_TurbSST final : public CUpwScalar { -private: - using Base = CUpwScalar; - using Base::nDim; - using Base::V_i; - using Base::V_j; - using Base::a0; - using Base::a1; - using Base::Flux; - using Base::Jacobian_i; - using Base::Jacobian_j; - using Base::ScalarVar_i; - using Base::ScalarVar_j; - using Base::idx; - using Base::bounded_scalar; - - /*! - * \brief Adds any extra variables to AD - */ - void ExtraADPreaccIn() override {} - - /*! - * \brief SST specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - Flux[0] = a0*V_i[idx.Density()]*ScalarVar_i[0] + a1*V_j[idx.Density()]*ScalarVar_j[0]; - Flux[1] = a0*V_i[idx.Density()]*ScalarVar_i[1] + a1*V_j[idx.Density()]*ScalarVar_j[1]; - - Jacobian_i[0][0] = a0; Jacobian_i[0][1] = 0.0; - Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = a0; - - Jacobian_j[0][0] = a1; Jacobian_j[0][1] = 0.0; - Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = a1; - } - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] config - Definition of the particular problem. - */ - CUpwSca_TurbSST(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) - : CUpwScalar(val_nDim, val_nVar, config) { bounded_scalar = config->GetBounded_Turb(); } -}; diff --git a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp b/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp deleted file mode 100644 index a4b2bbe264b..00000000000 --- a/SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp +++ /dev/null @@ -1,346 +0,0 @@ -/*! - * \file turb_diffusion.hpp - * \brief Declarations of numerics classes for discretization of - * viscous fluxes in turbulence problems. - * \author F. Palacios, T. Economon - * \version 8.5.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2026, 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 . - */ -#pragma once - -#include "../scalar/scalar_diffusion.hpp" - -/*! - * \class CAvgGrad_TurbSA - * \brief Class for computing viscous term using average of gradients (Spalart-Allmaras Turbulence model). - * \ingroup ViscDiscr - * \author A. Bueno. - */ -template -class CAvgGrad_TurbSA final : public CAvgGrad_Scalar { -private: - using Base = CAvgGrad_Scalar; - using Base::Laminar_Viscosity_i; - using Base::Laminar_Viscosity_j; - using Base::Density_i; - using Base::Density_j; - using Base::ScalarVar_i; - using Base::ScalarVar_j; - using Base::Proj_Mean_GradScalarVar; - using Base::proj_vector_ij; - using Base::Flux; - using Base::Jacobian_i; - using Base::Jacobian_j; - - const su2double sigma = 2.0/3.0; - const su2double cb2 = 0.622; - - const bool use_accurate_jacobians; - - /*! - * \brief Adds any extra variables to AD - */ - void ExtraADPreaccIn() override {} - - /*! - * \brief SA specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; - - /*--- Compute mean effective viscosity ---*/ - - /*--- First Term. Normal diffusion, and conservative part of the quadratic diffusion. - * ||grad nu_t||^2 = div(nu_t grad nu_t) - nu_t div grad nu_t ---*/ - const su2double nu_i = Laminar_Viscosity_i/Density_i; - const su2double nu_j = Laminar_Viscosity_j/Density_j; - const su2double nu_e = 0.5 * (nu_i + nu_j + (1 + cb2) * (ScalarVar_i[0] + ScalarVar_j[0])); - const su2double term_1 = nu_e; - - /* Second Term (quadratic diffusion, non conservative). */ - const su2double nu_tilde_i = ScalarVar_i[0]; - const su2double term_2 = cb2 * nu_tilde_i; - - const su2double diffusion_coefficient = term_1 - term_2; - Flux[0] = diffusion_coefficient * Proj_Mean_GradScalarVar[0] / sigma; - - if (implicit) { - /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/ - Jacobian_i[0][0] = -diffusion_coefficient * proj_vector_ij / sigma; - Jacobian_j[0][0] = diffusion_coefficient * proj_vector_ij / sigma; - - if (use_accurate_jacobians) { - /*--- The diffusion coefficient is also a function of nu_t. ---*/ - const su2double dTerm1_dnut_i = (1 + cb2) * 0.5; - const su2double dTerm1_dnut_j = (1 + cb2) * 0.5; - - const su2double dTerm2_dnut_i = cb2; - const su2double dTerm2_dnut_j = 0.0; - - const su2double dDC_dnut_i = dTerm1_dnut_i - dTerm2_dnut_i; - const su2double dDC_dnut_j = dTerm1_dnut_j - dTerm2_dnut_j; - - Jacobian_i[0][0] += dDC_dnut_i * Proj_Mean_GradScalarVar[0] / sigma; - Jacobian_j[0][0] += dDC_dnut_j * Proj_Mean_GradScalarVar[0] / sigma; - } - } - } - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] correct_grad - Whether to correct gradient for skewness. - * \param[in] config - Definition of the particular problem. - */ - CAvgGrad_TurbSA(unsigned short val_nDim, unsigned short val_nVar, - bool correct_grad, const CConfig* config) - : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config), - use_accurate_jacobians(config->GetUse_Accurate_Turb_Jacobians()) {} -}; - -/*! - * \class CAvgGrad_TurbSA_Neg - * \brief Class for computing viscous term using average of gradients (Spalart-Allmaras Turbulence model). - * \ingroup ViscDiscr - * \author F. Palacios - */ -template -class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar { -private: - using Base = CAvgGrad_Scalar; - using Base::Laminar_Viscosity_i; - using Base::Laminar_Viscosity_j; - using Base::Density_i; - using Base::Density_j; - using Base::ScalarVar_i; - using Base::ScalarVar_j; - using Base::Proj_Mean_GradScalarVar; - using Base::proj_vector_ij; - using Base::Flux; - using Base::Jacobian_i; - using Base::Jacobian_j; - - const su2double sigma = 2.0/3.0; - const su2double cn1 = 16.0; - const su2double cb2 = 0.622; - - /*! - * \brief Adds any extra variables to AD - */ - void ExtraADPreaccIn() override {} - - /*! - * \brief SA-neg specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; - - /*--- Compute mean effective viscosity ---*/ - - const su2double nu_i = Laminar_Viscosity_i/Density_i; - const su2double nu_j = Laminar_Viscosity_j/Density_j; - - const su2double nu_ij = 0.5 * (nu_i + nu_j); - const su2double nu_tilde_i = ScalarVar_i[0]; - const su2double nu_tilde_j = ScalarVar_j[0]; - const su2double nu_tilde_ij = 0.5 * (nu_tilde_i + nu_tilde_j); - - /*--- Following Diskin's implementation from 10.2514/1.J064629, they propose a new fn function - * to be evaluated at the cell to maintain positivity in the diffusion coefficient, which is - * used in both terms. The new fn term averaged across the face reverts to the original fn - * function. ---*/ - - /*--- Second Term (LHS) ---*/ - const su2double zeta_i = ((1 + cb2) * nu_tilde_ij - cb2 * nu_tilde_i) / nu_ij; - su2double fn_i = 1.0; - if (zeta_i < 0.0) { - fn_i = (cn1 + pow(zeta_i,3)) / (cn1 - pow(zeta_i,3)); - } - - const su2double term_1 = (nu_ij + (1 + cb2) * nu_tilde_ij * fn_i); - const su2double term_2 = cb2 * nu_tilde_i * fn_i; - Flux[0] = (term_1 - term_2) * Proj_Mean_GradScalarVar[0] / sigma; - - /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients - * Exact Jacobians were tested on multiple cases but resulted in divergence of all - * simulations, hence only frozen diffusion coefficient (approximate) Jacobians are used. ---*/ - - if (implicit) { - const su2double diffusion_coefficient = (term_1 - term_2); - - const su2double dGrad_dnut_i = -proj_vector_ij; - const su2double dGrad_dnut_j = proj_vector_ij; - - Jacobian_i[0][0] = diffusion_coefficient * dGrad_dnut_i / sigma; - Jacobian_j[0][0] = diffusion_coefficient * dGrad_dnut_j / sigma; - } - } - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] correct_grad - Whether to correct gradient for skewness. - * \param[in] config - Definition of the particular problem. - */ - CAvgGrad_TurbSA_Neg(unsigned short val_nDim, unsigned short val_nVar, - bool correct_grad, const CConfig* config) - : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config) {} -}; - -/*! - * \class CAvgGrad_TurbSST - * \brief Class for computing viscous term using average of gradient with correction (Menter SST turbulence model). - * \ingroup ViscDiscr - * \author A. Bueno. - */ -template -class CAvgGrad_TurbSST final : public CAvgGrad_Scalar { -private: - using Base = CAvgGrad_Scalar; - using Base::Laminar_Viscosity_i; - using Base::Laminar_Viscosity_j; - using Base::Eddy_Viscosity_i; - using Base::Eddy_Viscosity_j; - using Base::Density_i; - using Base::Density_j; - using Base::ScalarVar_i; - using Base::ScalarVar_j; - using Base::Proj_Mean_GradScalarVar; - using Base::proj_vector_ij; - using Base::Flux; - using Base::Jacobian_i; - using Base::Jacobian_j; - - const su2double sigma_k1; /*!< \brief Constants for the viscous terms, k-w (1), k-eps (2)*/ - const su2double sigma_k2; - const su2double sigma_om1; - const su2double sigma_om2; - const bool use_accurate_jacobians; - - su2double F1_i, F1_j; /*!< \brief Menter's first blending function */ - - /*! - * \brief Adds any extra variables to AD - */ - void ExtraADPreaccIn() override { - AD::SetPreaccIn(F1_i, F1_j); - } - - /*! - * \brief SST specific steps in the ComputeResidual method - * \param[in] config - Definition of the particular problem. - */ - void FinishResidualCalc(const CConfig* config) override { - const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; - - /*--- Compute the blended constant for the viscous terms ---*/ - const su2double sigma_kine_i = F1_i*sigma_k1 + (1.0 - F1_i)*sigma_k2; - const su2double sigma_kine_j = F1_j*sigma_k1 + (1.0 - F1_j)*sigma_k2; - const su2double sigma_omega_i = F1_i*sigma_om1 + (1.0 - F1_i)*sigma_om2; - const su2double sigma_omega_j = F1_j*sigma_om1 + (1.0 - F1_j)*sigma_om2; - - /*--- Compute mean effective dynamic viscosity ---*/ - const su2double diff_i_kine = Laminar_Viscosity_i + sigma_kine_i*Eddy_Viscosity_i; - const su2double diff_j_kine = Laminar_Viscosity_j + sigma_kine_j*Eddy_Viscosity_j; - const su2double diff_i_omega = Laminar_Viscosity_i + sigma_omega_i*Eddy_Viscosity_i; - const su2double diff_j_omega = Laminar_Viscosity_j + sigma_omega_j*Eddy_Viscosity_j; - - const su2double diff_kine = 0.5*(diff_i_kine + diff_j_kine); - const su2double diff_omega_T1 = 0.5*(diff_i_omega + diff_j_omega); - - /*--- We aim to treat the cross-diffusion as a diffusion term rather than a source term. - * Re-writing the cross-diffusion contribution as λ/w ∇w ∇k, where λ = (2 (1- F1) ρ σ_ω2) - * and expanding using the product rule for divergence theorem gives: ∇(w λ/w ∇k) - w ∇(λ/w ∇k). - * Discretising using FVM, gives: (λ)_ij ∇k - w_c (λ/w)_ij ∇k. where w_c is the cell centre value ---*/ - - const su2double lambda_i = 2 * (1 - F1_i) * Density_i * sigma_omega_i; - const su2double lambda_j = 2 * (1 - F1_j) * Density_j * sigma_omega_j; - const su2double lambda_ij = 0.5 * (lambda_i + lambda_j); - const su2double w_ij = 0.5 * (ScalarVar_i[1] + ScalarVar_j[1]); - - const su2double diff_omega_T2 = lambda_ij; - - const su2double diff_omega_T3 = -ScalarVar_i[1] * lambda_ij/w_ij; - - Flux[0] = diff_kine*Proj_Mean_GradScalarVar[0]; - Flux[1] = diff_omega_T1*Proj_Mean_GradScalarVar[1] + (diff_omega_T2 + diff_omega_T3)*Proj_Mean_GradScalarVar[0]; - - /*--- For Jacobians -> Use of TSL (Thin Shear Layer) approx. to compute derivatives of the gradients ---*/ - if (implicit) { - const su2double proj_on_rho_i = proj_vector_ij/Density_i; - const su2double proj_on_rho_j = proj_vector_ij/Density_j; - Jacobian_i[0][0] = -diff_kine*proj_on_rho_i; - Jacobian_i[0][1] = 0.0; - Jacobian_i[1][0] = (diff_omega_T2+diff_omega_T3)*-proj_on_rho_i; - Jacobian_i[1][1] = -diff_omega_T1*proj_on_rho_i; - - Jacobian_j[0][0] = diff_kine*proj_on_rho_j; - Jacobian_j[0][1] = 0.0; - Jacobian_j[1][0] = (diff_omega_T2+diff_omega_T3)*proj_on_rho_j; - Jacobian_j[1][1] = diff_omega_T1*proj_on_rho_j; - - if (use_accurate_jacobians) { - Jacobian_i[0][0] = -diff_kine*proj_on_rho_i; - Jacobian_i[0][1] = 0.0; - Jacobian_i[1][0] = (diff_omega_T2 + diff_omega_T3)*-proj_on_rho_i; - Jacobian_i[1][1] = -proj_on_rho_i * diff_omega_T1 - 2*lambda_ij*ScalarVar_j[1]/pow(ScalarVar_i[1]+ScalarVar_j[1],2) * Proj_Mean_GradScalarVar[0]; - - Jacobian_j[0][0] = diff_kine*proj_on_rho_j; - Jacobian_j[0][1] = 0.0; - Jacobian_j[1][0] = (diff_omega_T2 + diff_omega_T3)*proj_on_rho_j; - Jacobian_j[1][1] = proj_on_rho_j * diff_omega_T1 + 2*lambda_ij*ScalarVar_i[1]/pow(ScalarVar_i[1]+ScalarVar_j[1],2) * Proj_Mean_GradScalarVar[0]; - } - } - } - -public: - /*! - * \brief Constructor of the class. - * \param[in] val_nDim - Number of dimensions of the problem. - * \param[in] val_nVar - Number of variables of the problem. - * \param[in] constants - Constants of the model. - * \param[in] correct_grad - Whether to correct gradient for skewness. - * \param[in] config - Definition of the particular problem. - */ - CAvgGrad_TurbSST(unsigned short val_nDim, unsigned short val_nVar, - const su2double* constants, bool correct_grad, const CConfig* config) - : CAvgGrad_Scalar(val_nDim, val_nVar, correct_grad, config), - sigma_k1(constants[0]), - sigma_k2(constants[1]), - sigma_om1(constants[2]), - sigma_om2(constants[3]), - use_accurate_jacobians(config->GetUse_Accurate_Turb_Jacobians()) { - } - - /*! - * \brief Sets value of first blending function. - */ - void SetF1blending(su2double val_F1_i, su2double val_F1_j) override { - F1_i = val_F1_i; F1_j = val_F1_j; - } -}; diff --git a/SU2_CFD/include/numerics/turbulent/turb_sa_edge_flux.hpp b/SU2_CFD/include/numerics/turbulent/turb_sa_edge_flux.hpp new file mode 100644 index 00000000000..4e4b9594ce7 --- /dev/null +++ b/SU2_CFD/include/numerics/turbulent/turb_sa_edge_flux.hpp @@ -0,0 +1,150 @@ +/*! + * \file turb_sa_edge_flux.hpp + * \brief Spalart-Allmaras model as a third-layer scalar flux, see numerics/scalar/scalar_edge_flux.hpp. + * \author P. Gomes + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, 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 . + */ + +#pragma once + +#include "../scalar/scalar_edge_flux.hpp" + +/*! + * \class CScalarFlux_SA + * \ingroup ConvDiscr + * \ingroup ViscDiscr + * \brief Convection and diffusion of the Spalart-Allmaras model, non-conservative and with a + * diagonal (but asymmetric) diffusion coefficient. + * \note SA writes its own convective term rather than using the inherited CUpwScalarFlux one, + * because with stochastic backscatter active (nVar 4) the three Langevin equations are + * advected with a centered flux, unlike the plain upwind SA equation itself. + */ +template +class CScalarFlux_SA + : public CUpwScalarBase, FlowIndices, nDim, nVar> { + public: + static constexpr bool Conservative = false; + static constexpr bool DiagonalDiffusion = true; + static constexpr bool DiffusionReadsDensity = true; /*!< \brief The kinematic viscosities below. */ + + using Base = CUpwScalarBase; + using Int = typename Base::Int; + using Base::Base; + + private: + static constexpr passivedouble sigma = 2.0 / 3.0; /*!< \brief Constant of the diffusion term. */ + static constexpr passivedouble cb2 = 0.622; /*!< \brief Constant of the diffusion term. */ + + public: + /*! + * \brief SA convection, plus the centered advection of the backscatter equations when nVar > 1. + */ + template + FORCEINLINE void finalizeFlux(const FlowIndices&, const ScalarFluxOptions& opt, Int, const EdgeSide&, + Int, const EdgeSide&, const Double& a0, const Double& a1, + const CPair&, const CPair>& phi, + EdgeResidual& res) const { + const Double flux = a0 * phi.i.all(0) + a1 * phi.j.all(0); + + res.flux_i(0) += flux; + if (!opt.oneSided) res.flux_j(0) -= flux; + + if (opt.implicit) { + res.jac_ii(0, 0) += a0; + if (!opt.oneSided) { + res.jac_ij(0, 0) += a1; + res.jac_ji(0, 0) -= a0; + res.jac_jj(0, 0) -= a1; + } + } + + /*--- Stochastic backscatter: three Langevin equations, advected with the mean of the two + * upwinding weights and with no diffusion. ---*/ + const Double avg = 0.5 * (a0 + a1); + for (size_t iVar = 1; iVar < res.nVar; ++iVar) { + const Double flux_bs = avg * (phi.i.all(iVar) + phi.j.all(iVar)); + + res.flux_i(iVar) += flux_bs; + if (!opt.oneSided) res.flux_j(iVar) -= flux_bs; + + if (opt.implicit) { + res.jac_ii(iVar, iVar) += avg; + if (!opt.oneSided) { + res.jac_ij(iVar, iVar) += avg; + res.jac_ji(iVar, iVar) -= avg; + res.jac_jj(iVar, iVar) -= avg; + } + } + } + } + + /*! + * \brief Diffusion coefficients of both orientations of the edge. + * \note The coefficient is not symmetric: it uses the transported variable of the row it is + * going to be used for (the quadratic, non-conservative part of the diffusion term). + * Coefficients past index 0 are left at zero, the backscatter equations have no diffusion. + */ + template + FORCEINLINE CPair> coefficients(const FlowIndices& idx, Int iPoint, + const EdgeSide& side_i, Int jPoint, + const EdgeSide& side_j, + const CPair& rho) const { + const Double nu_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.LaminarViscosity()) / rho.i; + const Double nu_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.LaminarViscosity()) / rho.j; + + const Double nuTilde_i = gatherVariables(iPoint, side_i.scalarNodes.GetSolution(), 0); + const Double nuTilde_j = gatherVariables(jPoint, side_j.scalarNodes.GetSolution(), 0); + + const Double nu_e = 0.5 * (nu_i + nu_j + (1.0 + cb2) * (nuTilde_i + nuTilde_j)); + + Vector D_i, D_j; + D_i(0) = (nu_e - cb2 * nuTilde_i) / sigma; + D_j(0) = (nu_e - cb2 * nuTilde_j) / sigma; + for (size_t iVar = 1; iVar < nVar; ++iVar) { + D_i(iVar) = 0.0; + D_j(iVar) = 0.0; + } + return {D_i, D_j}; + } + + /*! + * \brief Extra Jacobian terms from the dependence of the diffusion coefficient on nu_tilde. + */ + template + FORCEINLINE void coefficientJacobians(const ScalarFluxOptions& opt, const Coefficients&, + const Vector& projGrad, EdgeResidual& res) const { + /*--- d(diffusion coefficient of i)/d(nu_tilde_i), and its counterpart w.r.t. nu_tilde_j; + * the coefficient of j is the same expression with i and j swapped, so the same two + * derivatives apply to both orientations. Both are per-edge constants, so the coefficients + * themselves are not read here. ---*/ + const Double dDC_dNuTilde_i = ((1.0 + cb2) * 0.5 - cb2) / sigma; + const Double dDC_dNuTilde_j = (1.0 + cb2) * 0.5 / sigma; + + res.jac_ii(0, 0) -= dDC_dNuTilde_i * projGrad(0); + if (opt.oneSided) return; + + res.jac_ij(0, 0) -= dDC_dNuTilde_j * projGrad(0); + res.jac_ji(0, 0) += dDC_dNuTilde_j * projGrad(0); + res.jac_jj(0, 0) += dDC_dNuTilde_i * projGrad(0); + } +}; diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 2c0cb51a8d2..9f72b6e37cc 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -1014,7 +1014,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { Residual[0] -= dk * Volume; Residual[1] -= dw * Volume; - /*--- Cross diffusion is included in the viscous fluxes, discretisation in turb_diffusion.hpp ---*/ + /*--- Cross diffusion is included in the viscous fluxes, not this source term. ---*/ /*--- Contribution due to 2D axisymmetric formulation ---*/ diff --git a/SU2_CFD/include/numerics/turbulent/turb_sst_edge_flux.hpp b/SU2_CFD/include/numerics/turbulent/turb_sst_edge_flux.hpp new file mode 100644 index 00000000000..b891b6e14d5 --- /dev/null +++ b/SU2_CFD/include/numerics/turbulent/turb_sst_edge_flux.hpp @@ -0,0 +1,147 @@ +/*! + * \file turb_sst_edge_flux.hpp + * \brief Menter SST model as a third-layer scalar flux, see numerics/scalar/scalar_edge_flux.hpp. + * \author P. Gomes + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, 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 . + */ + +#pragma once + +#include "../scalar/scalar_edge_flux.hpp" + +/*! + * \class CScalarFlux_SST + * \ingroup ViscDiscr + * \brief Convection and diffusion of the Menter SST model, conservative with a coupled (but + * neither symmetric nor diagonal) 2x2 diffusion matrix. + * \note SST writes no finalizeFlux of its own: the inherited CUpwScalarFlux one is exactly + * flux(iVar) = a0*rho_i*phi_i(iVar) + a1*rho_j*phi_j(iVar), Conservative weighting by + * density, which is the model's whole convective term. + */ +template +class CScalarFlux_SST + : public CUpwScalarBase, FlowIndices, nDim, nVar> { + public: + static constexpr bool Conservative = true; + static constexpr bool DiagonalDiffusion = false; + + using Base = CUpwScalarBase; + using Int = typename Base::Int; + using Base::Base; + + private: + /*--- Fixed regardless of SST_OPTIONS::version: only the production-limiter and source-term + * constants (alfa/gamma) differ by version, not these. ---*/ + static constexpr passivedouble sigma_k1 = 0.85; + static constexpr passivedouble sigma_k2 = 1.0; + static constexpr passivedouble sigma_om1 = 0.5; + static constexpr passivedouble sigma_om2 = 0.856; + + public: + /*! + * \brief Diffusion coefficients of both orientations of the edge, and the terms of the cross + * diffusion that the Jacobian correction below needs, so that neither the gathers nor + * the blending are repeated for it. + * \note The cross term reads the transported omega of whichever point its row is being written + * for, so it is not symmetric: i, read by i's row, uses omega at i; j, read by j's row, + * uses omega at j. Every other entry is an i/j average, so it is the same in both. + */ + struct CCoefficients { + Matrix i, j; + Double lambda_ij, omega_i, omega_j; + }; + + template + FORCEINLINE CCoefficients coefficients(const FlowIndices& idx, Int iPoint, const EdgeSide& side_i, + Int jPoint, const EdgeSide& side_j, + const CPair& rho) const { + const Double mu_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.LaminarViscosity()); + const Double mu_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.LaminarViscosity()); + const Double muT_i = gatherVariables(iPoint, side_i.flowNodes->GetPrimitive(), idx.EddyViscosity()); + const Double muT_j = gatherVariables(jPoint, side_j.flowNodes->GetPrimitive(), idx.EddyViscosity()); + + const Double F1_i = gatherVariables(iPoint, side_i.scalarNodes.GetF1blending()); + const Double F1_j = gatherVariables(jPoint, side_j.scalarNodes.GetF1blending()); + + CCoefficients D; + D.omega_i = gatherVariables(iPoint, side_i.scalarNodes.GetSolution(), 1); + D.omega_j = gatherVariables(jPoint, side_j.scalarNodes.GetSolution(), 1); + + const Double sigma_kine_i = F1_i * sigma_k1 + (1.0 - F1_i) * sigma_k2; + const Double sigma_kine_j = F1_j * sigma_k1 + (1.0 - F1_j) * sigma_k2; + const Double sigma_omega_i = F1_i * sigma_om1 + (1.0 - F1_i) * sigma_om2; + const Double sigma_omega_j = F1_j * sigma_om1 + (1.0 - F1_j) * sigma_om2; + + const Double diff_kine = 0.5 * ((mu_i + sigma_kine_i * muT_i) + (mu_j + sigma_kine_j * muT_j)); + const Double diff_omega = 0.5 * ((mu_i + sigma_omega_i * muT_i) + (mu_j + sigma_omega_j * muT_j)); + + const Double lambda_i = 2.0 * (1.0 - F1_i) * rho.i * sigma_omega_i; + const Double lambda_j = 2.0 * (1.0 - F1_j) * rho.j * sigma_omega_j; + D.lambda_ij = 0.5 * (lambda_i + lambda_j); + const Double w_ij = 0.5 * (D.omega_i + D.omega_j); + + /*--- Cross-diffusion coefficient: a divergence-theorem term (diff_omega_T2) plus a cell + * centre correction (diff_omega_T3) that reads the transported omega of the row's own point. ---*/ + const Double diff_omega_T2 = D.lambda_ij; + const Double diff_omega_T3_i = -D.omega_i * D.lambda_ij / w_ij; + const Double diff_omega_T3_j = -D.omega_j * D.lambda_ij / w_ij; + + /*--- D.i(0,1) and D.j(0,1) are left zero: there is no diffusive coupling from omega into + * the k row. ---*/ + D.i = Double(0.0); + D.j = Double(0.0); + D.i(0, 0) = diff_kine; + D.i(1, 1) = diff_omega; + D.i(1, 0) = diff_omega_T2 + diff_omega_T3_i; + + D.j(0, 0) = diff_kine; + D.j(1, 1) = diff_omega; + D.j(1, 0) = diff_omega_T2 + diff_omega_T3_j; + + return D; + } + + /*! + * \brief Extra Jacobian terms from the dependence of the cross-diffusion coefficient on omega. + * \note diff_omega_T3_i and diff_omega_T3_j both depend on omega_i and omega_j through w_ij, so + * each of the four blocks needs a correction beyond the one diffusionTerms already applies + * through projGrad. The correction only depends on which point's omega is being + * differentiated against, not on which row it lands in: differentiating against omega_i + * gives +E_j in both jac_ii and jac_ji, differentiating against omega_j gives -E_i in both + * jac_ij and jac_jj. + */ + template + FORCEINLINE void coefficientJacobians(const ScalarFluxOptions& opt, const CCoefficients& D, + const Vector& projGrad, EdgeResidual& res) const { + const Double denom = pow(D.omega_i + D.omega_j, 2.0); + const Double E_i = 2.0 * D.lambda_ij * D.omega_i / denom * projGrad(0); + const Double E_j = 2.0 * D.lambda_ij * D.omega_j / denom * projGrad(0); + + res.jac_ii(1, 1) += E_j; + if (opt.oneSided) return; + + res.jac_ij(1, 1) -= E_i; + res.jac_ji(1, 1) += E_j; + res.jac_jj(1, 1) -= E_i; + } +}; diff --git a/SU2_CFD/include/numerics/util.hpp b/SU2_CFD/include/numerics/util.hpp new file mode 100644 index 00000000000..92e5144b9c9 --- /dev/null +++ b/SU2_CFD/include/numerics/util.hpp @@ -0,0 +1,638 @@ +/*! + * \file util.hpp + * \brief Generic auxiliary functions. + * \author P. Gomes + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, 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 . + */ + +#pragma once + +#include + +#include "../../../Common/include/option_structure.hpp" +#include "../../../Common/include/parallelization/vectorization.hpp" +#include "../../../Common/include/containers/C2DContainer.hpp" +#include "../../../Common/include/linear_algebra/CSysVector.hpp" +#include "../../../Common/include/linear_algebra/CSysMatrix.hpp" + +/*! + * \enum UpdateType + * \brief Ways to update vectors and system matrices. + * COLORING is the typical i/j update, whereas for REDUCTION + * the fluxes are stored and the matrix diagonal is not modified. + */ +enum class UpdateType { COLORING, REDUCTION }; + +#ifdef CODI_FORWARD_TYPE +using SparseMatrixType = CSysMatrix; +#else +using SparseMatrixType = CSysMatrix; +#endif + +/*! + * \brief Alignment of the static containers backing a flux value type. + * \note Yields the type's own alignment for a SIMD array, and the plain type's natural + * alignment for a scalar; C2DContainer's AlignSize also accepts 0 to mean its own + * default, but that reaches `alignas(0)` on the static specializations, which some + * compilers warn about even though it is a no-op, so a real value is passed instead. + */ +template +struct CAlignTraits { + enum : size_t { Align = alignof(Type) }; +}; + +template +struct CAlignTraits> { + enum : size_t { Align = simd::Array::Align }; +}; + +/*! + * \brief Static vector and matrix types. + * \note These should be used instead of C-style arrays. + */ +template +using Vector = C2DContainer::Align, Size, 1>; + +template +using Matrix = C2DContainer::Align, Rows, Cols>; + +/*! + * \brief The flux value type and lane count that go with an index type. + * \note There is exactly one floating type in play, su2double, active under AD; a plain + * integral index reads one of it, a lane-vector index reads a lane-vector of it. This + * is what lets every helper below deduce its value type from the index it is handed, + * instead of a caller naming it explicitly. + */ +template +struct CValueTraits; + +template <> +struct CValueTraits { + using Double = su2double; + static constexpr size_t Size = 1; +}; + +template +struct CValueTraits> { + using Double = simd::Array; + static constexpr size_t Size = N; +}; + +/*! + * \brief Index type and lane count that go with a flux value type, the converse of + * CValueTraits, used where the value type is already known (e.g. a class template + * parameter) and the index type is what needs deriving. + */ +template +struct CLaneTraits; + +template <> +struct CLaneTraits { + using Int = unsigned long; + static constexpr bool IsArray = false; +}; + +template +struct CLaneTraits> { + using Int = simd::Array; + static constexpr bool IsArray = true; +}; + +/*! + * \brief Constexpr version of max. + */ +inline constexpr size_t Max(size_t a, size_t b) { return a > b ? a : b; } + +/*! + * \brief Simple pair type for i/j variables. + */ +template +struct CPair { + T i, j; +}; + +/*! + * \brief Blocks a template parameter from participating in argument deduction. + * \note Deduction never applies a user conversion, so a parameter typed plain Double would + * force a caller passing a bare su2double constant (kappa, a limiter ramp) to have + * already broadcast it. Wrapping the parameter type here defers Double entirely to + * the other, genuinely deduced arguments, and the broadcast then happens as an + * ordinary implicit conversion at the call. + */ +template +struct CIdentity { + using type = T; +}; +template +using CNonDeduced = typename CIdentity::type; + +/*! + * \brief Equation count of a model whose value is only known at runtime. + */ +constexpr size_t Dynamic = size_t(-1); + +/*! + * \brief Backing size of the static arrays of a dynamic model. + * \note The scalar numerics cap the equation count at this value and error above it, so a + * configuration that fits them fits these kernels. + */ +constexpr size_t MaxScalarVar = 8; + +/*! + * \brief Residual of one edge, accumulated by the convective and the diffusive terms. + * \note flux_i and flux_j are the contributions to the rows of i and j. They are opposite + * for a conservative term and independent for a non-conservative one. The Jacobians + * map onto the ii, ij, ji and jj blocks of the edge. A dynamic model sizes the storage + * with the maximum and iterates to nVar, which the scheme sets from the solver. + */ +template +struct EdgeResidual { + /*!< \brief The Matrix a static nVar==1 model would otherwise need degenerates + * to vector-only indexing in C2DContainer (its RowMajor, one-row specialization), so the + * backing is never smaller than 2; the unused padding row/column is simply never visited, + * every loop here and in the model bounds itself to nVar, not Size. */ + static constexpr size_t Size = Max(2, (nVar_ == Dynamic) ? MaxScalarVar : nVar_); + + Vector flux_i, flux_j; + Matrix jac_ii, jac_ij, jac_ji, jac_jj; + const size_t nVar; + + /*! + * \brief Zero the terms of the equations in use, so that both terms can accumulate into them. + * \note A static model zeroes its whole storage with constant trip counts; a dynamic one + * zeroes the leading nVar rows and columns and leaves the rest of the backing untouched. + */ + FORCEINLINE explicit EdgeResidual(size_t nEqn) : nVar(nEqn) { + for (size_t iVar = 0; iVar < nVar; ++iVar) { + flux_i(iVar) = 0.0; + flux_j(iVar) = 0.0; + for (size_t jVar = 0; jVar < nVar; ++jVar) { + jac_ii(iVar, jVar) = 0.0; + jac_ij(iVar, jVar) = 0.0; + jac_ji(iVar, jVar) = 0.0; + jac_jj(iVar, jVar) = 0.0; + } + } + } +}; + +/*! + * \brief Dot product. + */ +template +FORCEINLINE auto dot(ForwardIterator iterator, const T* ptr) -> typename std::decay::type { + typename std::decay::type sum = 0.0; + for (size_t iDim = 0; iDim < nDim; ++iDim) { + sum += *(iterator++) * ptr[iDim]; + } + return sum; +} + +/*! + * \overload Dot product. + */ +template +FORCEINLINE Double dot(ForwardIterator iterator, const Vector& vector) { + return dot(iterator, vector.data()); +} + +/*! + * \overload Dot product. + */ +template +FORCEINLINE Double dot(const Vector& a, const Vector& b) { + return dot(a.data(), b.data()); +} + +/*! + * \brief Squared norm. + */ +template +FORCEINLINE auto squaredNorm(ForwardIterator iterator) -> typename std::decay::type { + typename std::decay::type sum = 0.0; + for (size_t iDim = 0; iDim < nDim; ++iDim) { + sum += pow(*(iterator++), 2); + } + return sum; +} + +/*! + * \overload Squared norm. + */ +template +FORCEINLINE Double squaredNorm(const Vector& vector) { + return squaredNorm(vector.data()); +} + +/*! + * \brief Tangential projection. + */ +template +FORCEINLINE Vector tangentProjection(const Matrix& tensor, + const Vector& unitVector) { + Vector proj; + for (size_t iDim = 0; iDim < nDim; ++iDim) proj(iDim) = dot(tensor[iDim], unitVector); + + Double normalProj = dot(proj, unitVector); + + for (size_t iDim = 0; iDim < nDim; ++iDim) proj(iDim) -= normalProj * unitVector(iDim); + + return proj; +} + +/*! + * \brief Vector norm. + */ +template +FORCEINLINE Double norm(const Vector& vector) { + return sqrt(squaredNorm(vector)); +} + +#ifndef CODI_REVERSE_TYPE +/*! + * \brief Gather a single variable, from column iVar (0 by default) of row iPoint of a + * 2D container, or from index iPoint of a 1D container. + */ +template ::Double> +FORCEINLINE Double gatherVariables(Int iPoint, const Container& vars, size_t iVar = 0) { + return vars.template get>(iPoint, iVar)(0); +} + +/*! + * \brief Gather nVar contiguous variables starting at column iVar (0 by default) of row + * iPoint of a 2D container. + */ +template ::Double> +FORCEINLINE Vector gatherVariables(Int iPoint, const Container& vars, size_t iVar = 0) { + return vars.template get>(iPoint, iVar); +} + +/*! + * \brief Gather an nRows x nCols block of a 3D container, from outer index iPoint and + * starting at middle index iRow. + */ +template ::Double> +FORCEINLINE Matrix gatherVariables(Int iPoint, const Container& vars, size_t iRow = 0) { + return vars.template get>(iPoint, iRow); +} +#else + +namespace { +/*--- Rank of a container, from the element accessors it offers: a 1D container is indexed by + * point alone, a 3D one by (point, row, column). A 3D container also answers a two-argument + * call, but with an offset sub-matrix view rather than a scalar, which is why the rank is + * detected up front instead of trying the access forms in turn. ---*/ +template +struct Is1D : std::false_type {}; +template +struct Is1D()(0ul))>> : std::true_type {}; + +template +struct Is3D : std::false_type {}; +template +struct Is3D()(0ul, 0ul, 0ul))>> : std::true_type {}; + +/*--- One lane of an index or of a gathered value, the whole thing when there are no lanes. ---*/ +FORCEINLINE unsigned long lane(unsigned long iPoint, size_t) { return iPoint; } +template +FORCEINLINE unsigned long lane(const simd::Array& iPoint, size_t k) { return iPoint[k]; } +FORCEINLINE su2double& lane(su2double& x, size_t) { return x; } +template +FORCEINLINE T& lane(simd::Array& x, size_t k) { return x[k]; } + +/*--- Register one source element as a preaccumulation input, passing it through for the copy. + * The registration has to happen here, on the reference into the container's own storage: a + * copy has a fresh identifier of its own, and registering that instead would sever the source + * from the statement EndPreacc() stores. ---*/ +FORCEINLINE const su2double& preaccIn(const su2double& value) { + AD::SetPreaccIn(value); + return value; +} +} // namespace + +template ::Double> +FORCEINLINE Double gatherVariables(Int iPoint, const Container& vars, size_t iVar = 0) { + Double x; + for (size_t k = 0; k < CValueTraits::Size; ++k) { + if constexpr (Is1D::value) + lane(x, k) = preaccIn(vars(lane(iPoint, k))); + else + lane(x, k) = preaccIn(vars(lane(iPoint, k), iVar)); + } + return x; +} + +template ::Double> +FORCEINLINE Vector gatherVariables(Int iPoint, const Container& vars, size_t iVar = 0) { + Vector x; + for (size_t i = 0; i < nVar; ++i) { + for (size_t k = 0; k < CValueTraits::Size; ++k) { + if constexpr (Is3D::value) + lane(x(i), k) = preaccIn(vars(lane(iPoint, k), iVar, i)); + else + lane(x(i), k) = preaccIn(vars(lane(iPoint, k), iVar + i)); + } + } + return x; +} + +template ::Double> +FORCEINLINE Matrix gatherVariables(Int iPoint, const Container& vars, size_t iRow = 0) { + Matrix x; + for (size_t i = 0; i < nRows; ++i) { + for (size_t j = 0; j < nCols; ++j) { + for (size_t k = 0; k < CValueTraits::Size; ++k) { + lane(x(i, j), k) = preaccIn(vars(lane(iPoint, k), iRow + i, j)); + } + } + } + return x; +} +#endif + +/*! + * \brief Register the leading nVar entries of a static vector as preaccumulation outputs. + * \note A lane vector is registered one lane at a time, a scalar in one call; a kernel therefore + * reaches this rather than AD::SetPreaccOut directly, and reads the same whichever value + * type it is bound to. + */ +template +FORCEINLINE void setPreaccOut(Vector& x, size_t nVar) { + if constexpr (!CLaneTraits::IsArray) { + AD::SetPreaccOut(x, static_cast(nVar)); + } else { + AD::SetPreaccOut(x, static_cast(nVar), Double::Size); + } +} + +/*! + * \brief Stop the AD preaccumulation. + */ +template +FORCEINLINE void stopPreacc(Vector& x) { + setPreaccOut(x, nVar); + AD::EndPreacc(); +} + +/*! + * \brief Distance vector, from point i to point j of one container. + */ +template ::Double> +FORCEINLINE Vector distanceVector(Int iPoint, Int jPoint, const Container& coords) { + return distanceVector(iPoint, coords, jPoint, coords); +} + +/*! + * \brief Distance vector, from point i of one container to point j of another. + * \note The two endpoints of a boundary flux read different containers, the solver's own + * and the marker's ghost one; the interior edge loop passes the same container twice. + */ +template ::Double> +FORCEINLINE Vector distanceVector(Int iPoint, const Container& coords_i, Int jPoint, + const Container& coords_j) { + auto coord_i = gatherVariables(iPoint, coords_i); + auto coord_j = gatherVariables(jPoint, coords_j); + Vector vector_ij; + for (size_t iDim = 0; iDim < nDim; ++iDim) { + vector_ij(iDim) = coord_j(iDim) - coord_i(iDim); + } + return vector_ij; +} + +/*! + * \brief Blended difference for U-MUSCL reconstruction. + * \param[in] gradProj - Gradient projection at point i: dot(grad_i, vector_ij). + * \param[in] delta - Centered difference: V_j - V_i. + * \param[in] kappa - Blending parameter. + * \return Blended difference for reconstruction from point i. + */ +template +FORCEINLINE Double umusclProjection(const Double& gradProj, const Double& delta, const CNonDeduced& kappa) { + /*-------------------------------------------------------------------*/ + /*--- The MUSCL kappa-scheme reconstruction is typically written: ---*/ + /*--- V_L = V_i + 0.25 * dV_ij^kap, where ---*/ + /*--- dV_ij^kap = (1-kappa) dV_ij^upw + (1+kappa) dV_ij^cen, ---*/ + /*--- dV_ij^cen = V_j - V_i, ---*/ + /*--- dV_ij^upw = 2 grad(Vi) dot vector_ij - dV_ij^cen. ---*/ + /*--- To maintain proper scaling for edge limiters, the result of ---*/ + /*--- this function is 0.5 * dV_ij^kap. ---*/ + /*-------------------------------------------------------------------*/ + return (1.0 - kappa) * gradProj + kappa * delta; +} + +/*! + * \brief Reads the gradient rows of one point as one block, gathered up front. + * \note This is what a kernel whose variable count is a compile-time constant above one wants: + * one gather of the whole nVarGrad x nDim block instead of nVarGrad of them. + */ +template +struct CGradientBlock { + using Int = typename CLaneTraits::Int; + Matrix rows; + + FORCEINLINE CGradientBlock(Int iPoint, const Gradient_t& gradient, size_t iRow) + : rows(gatherVariables(iPoint, gradient, iRow)) {} + + FORCEINLINE Double project(size_t iVar, const Vector& vector_ij) const { + return dot(rows[iVar], vector_ij); + } +}; + +/*! + * \brief Reads the gradient rows of one point one row at a time. + * \note This is what a runtime variable count forces, since the block shape would not be a + * compile-time constant, and what a single variable forces, since its block would be a + * Matrix: that shape satisfies IsVector and would silently degenerate into + * a lone scalar instead of failing to compile. + */ +template +struct CGradientRows { + using Int = typename CLaneTraits::Int; + const Int iPoint; + const Gradient_t& gradient; + const size_t iRow; + + FORCEINLINE Double project(size_t iVar, const Vector& vector_ij) const { + return dot(gatherVariables(iPoint, gradient, iRow + iVar), vector_ij); + } +}; + +/*! + * \brief Gradient reader of one point, blocked or row by row according to nVarGrad_. + */ +template +FORCEINLINE auto gradientReader(typename CLaneTraits::Int iPoint, const Gradient_t& gradient, size_t iRow) { + if constexpr (nVarGrad_ > 1) { + return CGradientBlock(iPoint, gradient, iRow); + } else { + return CGradientRows{iPoint, gradient, iRow}; + } +} + +/*! + * \brief How the reconstructed differences are limited. + */ +enum class MusclLimiter { NONE, EDGE, POINT }; + +/*! + * \brief U-MUSCL reconstruction of nVarGrad variables, from the gradient rows starting at iRow. + * \note The limiter kind is a template parameter so that the choice is made once, outside the + * loop, by the dispatching overload below. + * \param[in] iRow - Starting row of gradient (and column of limiter) to read, for reconstructing + * a slice of a larger set of gradients (e.g. only the velocity out of the primitives). + * \param[in] nVarGradRuntime - Variable count of a Dynamic model, known only at runtime; ignored + * (falling back to nVarGrad_ or VarType::nVar) when left at its default of 0. + */ +template +FORCEINLINE void muscl(typename CLaneTraits::Int iPoint, typename CLaneTraits::Int jPoint, + const Vector& vector_ij, const Gradient_t& gradient, const Limiter_t& limiter, + size_t iRow, CPair& V, const CNonDeduced& kappa, + const CNonDeduced& umusclRamp, size_t nVarGradRuntime) { + const size_t nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : (nVarGradRuntime > 0 ? nVarGradRuntime : VarType::nVar); + + const auto grad_i = gradientReader(iPoint, gradient, iRow); + const auto grad_j = gradientReader(jPoint, gradient, iRow); + + for (size_t iVar = 0; iVar < nVarGrad; ++iVar) { + /*--- Centered difference, needed for the U-MUSCL projection and the edge limiter. ---*/ + const Double delta_ij = V.j.all(iVar) - V.i.all(iVar); + + /*--- U-MUSCL reconstructed differences, to be halved when applied. ---*/ + const Double proj_i = umusclRamp * umusclProjection(grad_i.project(iVar, vector_ij), delta_ij, kappa); + const Double proj_j = umusclRamp * umusclProjection(grad_j.project(iVar, vector_ij), delta_ij, kappa); + + Double lim_i = 1.0, lim_j = 1.0; + if constexpr (limiterKind == MusclLimiter::EDGE) { + const Double delta_ij_2 = pow(delta_ij, 2) + 1e-6; + lim_i = (delta_ij_2 + proj_i * delta_ij) / (pow(proj_i, 2) + delta_ij_2); + lim_j = (delta_ij_2 + proj_j * delta_ij) / (pow(proj_j, 2) + delta_ij_2); + } else if constexpr (limiterKind == MusclLimiter::POINT) { + lim_i = gatherVariables(iPoint, limiter, iRow + iVar); + lim_j = gatherVariables(jPoint, limiter, iRow + iVar); + } + + /*--- Apply reconstruction: V_L = V_i + 0.5 * lim * dV_ij^kap ---*/ + V.i.all(iVar) += 0.5 * lim_i * proj_i; + V.j.all(iVar) -= 0.5 * lim_j * proj_j; + } +} + +/*! + * \brief Reconstruct a slice of nVarGrad variables starting at column iRow, dispatching on the + * limiter type; shared by the flow and the scalar reconstructions so both call one body. + */ +template +FORCEINLINE void reconstruct(typename CLaneTraits::Int iPoint, typename CLaneTraits::Int jPoint, + const Vector& vector_ij, const Gradient_t& gradient, + const Limiter_t& limiter, LIMITER limiterType, size_t iRow, CPair& V, + const CNonDeduced& kappa, const CNonDeduced& umusclRamp, + size_t nVarGradRuntime = 0) { + switch (limiterType) { + case LIMITER::NONE: + muscl(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp, + nVarGradRuntime); + break; + case LIMITER::VAN_ALBADA_EDGE: + muscl(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp, + nVarGradRuntime); + break; + default: + muscl(iPoint, jPoint, vector_ij, gradient, limiter, iRow, V, kappa, umusclRamp, + nVarGradRuntime); + break; + } +} + +/*! + * \brief Update the matrix and right-hand-side of a linear system with one conservative flux. + */ +template ::Int> +FORCEINLINE void updateLinearSystem(Int iEdge, Int iPoint, Int jPoint, bool implicit, UpdateType updateType, + Double updateMask, const Vector& flux, + const Matrix& jac_i, const Matrix& jac_j, + CSysVector& vector, SparseMatrixType& matrix) { + if (updateType == UpdateType::COLORING) { + vector.UpdateBlocks(iPoint, jPoint, flux, updateMask); + if (implicit) { + auto wasActive = AD::BeginPassive(); + matrix.SetBlocks(iEdge, iPoint, jPoint, jac_i, jac_j, updateMask); + AD::EndPassive(wasActive); + } + } else { + vector.SetBlock(iEdge, flux, updateMask); + if (implicit) { + auto wasActive = AD::BeginPassive(); + matrix.SetBlocks(iEdge, jac_i, jac_j, updateMask); + AD::EndPassive(wasActive); + } + } +} + +/*! + * \brief Update the matrix and right-hand-side of a linear system with two independent row + * contributions and four independent Jacobian blocks. + * \note It carries a second CSysVector, the target of flux_j under UpdateType::REDUCTION and + * unused under COLORING, where both rows are written directly. + * \note The residual is the same under both update types, the Jacobian is not: REDUCTION writes + * the off-diagonal blocks only and CSysMatrix::SetDiagonalAsColumnSum then derives each + * diagonal block as minus the sum of its column, which equals the jac_ii and jac_jj computed + * here only where the flux is conservative (jac_ii == -jac_ji). A model whose diffusion + * coefficients differ between the two orientations of an edge, i.e. one that evaluates a + * non-conservative term at the point whose row it is writing, therefore converges along a + * slightly different path under the reducer, to the same solution. + */ +template ::Int> +FORCEINLINE void updateLinearSystem(Int iEdge, Int iPoint, Int jPoint, bool implicit, UpdateType updateType, + Double updateMask, const EdgeResidual& res, + CSysVector& vector, CSysVector& vectorDiff, + SparseMatrixType& matrix) { + if (updateType == UpdateType::COLORING) { + vector.AddBlock(iPoint, res.flux_i, updateMask); + vector.AddBlock(jPoint, res.flux_j, updateMask); + if (implicit) { + auto wasActive = AD::BeginPassive(); + matrix.SetBlocks(iEdge, iPoint, jPoint, res.jac_ii, res.jac_ij, res.jac_ji, res.jac_jj, updateMask); + AD::EndPassive(wasActive); + } + } else { + vector.SetBlock(iEdge, res.flux_i, updateMask); + vectorDiff.SetBlock(iEdge, res.flux_j, updateMask); + if (implicit) { + auto wasActive = AD::BeginPassive(); + matrix.SetOffDiagBlocks(iEdge, res.jac_ij, res.jac_ji, updateMask); + AD::EndPassive(wasActive); + } + } +} + +/*! + * \brief Store the (scalar) mass flux of an edge, e.g. for "bounded scalar" transport equations. + * \note No-op if "target" is null. As with CEdge's Nodes/Normal, edges within a SIMD group are + * contiguous (coloring groups are multiples of the SIMD size), so this is a plain vectorized store + * starting at iEdge[0], relying on "target" being padded to a multiple of the SIMD size. + */ +template ::Int> +FORCEINLINE void updateEdgeMassFlux(Int iEdge, const Double& massFlux, su2activevector* target) { + if (target) massFlux.store(&(*target)[iEdge[0]]); +} diff --git a/SU2_CFD/include/numerics_simd/CNumericsSIMD.hpp b/SU2_CFD/include/numerics_simd/CNumericsSIMD.hpp index a7726a6436a..076c08517b4 100644 --- a/SU2_CFD/include/numerics_simd/CNumericsSIMD.hpp +++ b/SU2_CFD/include/numerics_simd/CNumericsSIMD.hpp @@ -28,15 +28,7 @@ #pragma once #include "../../../Common/include/parallelization/vectorization.hpp" -#include "../../../Common/include/containers/C2DContainer.hpp" - -/*! - * \enum UpdateType - * \brief Ways to update vectors and system matrices. - * COLORING is the typical i/j update, whereas for REDUCTION - * the fluxes are stored and the matrix diagonal is not modified. - */ -enum class UpdateType {COLORING, REDUCTION}; +#include "../numerics/util.hpp" /*! * \brief Define Double and Int SIMD types. @@ -45,18 +37,10 @@ using Double = simd::Array; using Int = simd::Array; /*--- Forward declare a few classes used in name only by the interface. ---*/ -template class CSysVector; -template class CSysMatrix; class CConfig; class CGeometry; class CVariable; -#ifdef CODI_FORWARD_TYPE -using SparseMatrixType = CSysMatrix; -#else -using SparseMatrixType = CSysMatrix; -#endif - /*! * \class CNumericsSIMD * \ingroup ConvDiscr diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index f3669d45423..800167d4b0f 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -32,149 +32,6 @@ #include "../variables.hpp" #include "../../../variables/CNSVariable.hpp" -/*! - * \brief Blended difference for U-MUSCL reconstruction. - * \param[in] gradProj - Gradient projection at point i: dot(grad_i, vector_ij). - * \param[in] delta - Centered difference: V_j - V_i. - * \param[in] kappa - Blending parameter. - * \return Blended difference for reconstruction from point i. - */ -FORCEINLINE Double umusclProjection(const Double& gradProj, - const Double& delta, - const Double& kappa) { - /*-------------------------------------------------------------------*/ - /*--- The MUSCL kappa-scheme reconstruction is typically written: ---*/ - /*--- V_L = V_i + 0.25 * dV_ij^kap, where ---*/ - /*--- dV_ij^kap = (1-kappa) dV_ij^upw + (1+kappa) dV_ij^cen, ---*/ - /*--- dV_ij^cen = V_j - V_i, ---*/ - /*--- dV_ij^upw = 2 grad(Vi) dot vector_ij - dV_ij^cen. ---*/ - /*--- To maintain proper scaling for edge limiters, the result of ---*/ - /*--- this function is 0.5 * dV_ij^kap. ---*/ - /*-------------------------------------------------------------------*/ - return (1.0 - kappa) * gradProj + kappa * delta; -} - -/*! - * \brief MUSCL reconstruction of the specified variable. - * \note The result should be halved when added to i (or subtracted from j). - * \param[in] grad_i - Gradient vector at point i. - * \param[in] vector_ij - Distance vector from i to j. - * \param[in] delta - Centered difference: V_j - V_i. - * \param[in] iVar - Variable index. - * \param[in] kappa - Blending coefficient. - * \param[in] umusclRamp - MUSCL 1st-2nd order ramp times Newton-Krylov relaxation. - * \return Variable reconstructed from point i. - */ -template -FORCEINLINE Double musclReconstruction(const GradType& grad, - const VectorDbl& vector_ij, - const Double& delta, - const size_t iVar, - const Double& kappa, - const Double& umusclRamp) { - const Double proj = dot(grad[iVar], vector_ij); - return umusclRamp * umusclProjection(proj, delta, kappa); -} - -/*! - * \brief Unlimited reconstruction. - */ -template -FORCEINLINE void musclUnlimited(const Int& iPoint, - const Int& jPoint, - const VectorDbl& vector_ij, - const Gradient_t& gradient, - CPair& V, - const Double& kappa, - const Double& umusclRamp) { - constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar; - - auto grad_i = gatherVariables(iPoint, gradient); - auto grad_j = gatherVariables(jPoint, gradient); - - for (size_t iVar = 0; iVar < nVarGrad; ++iVar) { - /*--- Centered difference, needed for U-MUSCL projection ---*/ - const Double delta_ij = V.j.all(iVar) - V.i.all(iVar); - - /*--- U-MUSCL reconstructed variables ---*/ - const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, umusclRamp); - const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, umusclRamp); - - /*--- Apply reconstruction: V_L = V_i + 0.5 * dV_ij^kap ---*/ - V.i.all(iVar) += 0.5 * proj_i; - V.j.all(iVar) -= 0.5 * proj_j; - } -} - -/*! - * \brief Limited reconstruction with point-based limiter. - */ -template -FORCEINLINE void musclPointLimited(const Int& iPoint, - const Int& jPoint, - const VectorDbl& vector_ij, - const Limiter_t& limiter, - const Gradient_t& gradient, - CPair& V, - const Double& kappa, - const Double& umusclRamp) { - constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar; - - auto lim_i = gatherVariables(iPoint, limiter); - auto lim_j = gatherVariables(jPoint, limiter); - - auto grad_i = gatherVariables(iPoint, gradient); - auto grad_j = gatherVariables(jPoint, gradient); - - for (size_t iVar = 0; iVar < nVarGrad; ++iVar) { - /*--- Centered difference, needed for U-MUSCL projection ---*/ - const Double delta_ij = V.j.all(iVar) - V.i.all(iVar); - - /*--- U-MUSCL reconstructed variables ---*/ - const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, umusclRamp); - const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, umusclRamp); - - /*--- Apply reconstruction: V_L = V_i + 0.5 * lim * dV_ij^kap ---*/ - V.i.all(iVar) += 0.5 * lim_i(iVar) * proj_i; - V.j.all(iVar) -= 0.5 * lim_j(iVar) * proj_j; - } -} - -/*! - * \brief Limited reconstruction with edge-based limiter. - */ -template -FORCEINLINE void musclEdgeLimited(const Int& iPoint, - const Int& jPoint, - const VectorDbl& vector_ij, - const Gradient_t& gradient, - CPair& V, - const Double& kappa, - const Double& umusclRamp) { - constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar; - - auto grad_i = gatherVariables(iPoint, gradient); - auto grad_j = gatherVariables(jPoint, gradient); - - for (size_t iVar = 0; iVar < nVarGrad; ++iVar) { - /*--- Centered difference, needed for U-MUSCL projection and limiter ---*/ - const Double delta_ij = V.j.all(iVar) - V.i.all(iVar); - const Double delta_ij_2 = pow(delta_ij, 2) + 1e-6; - - /*--- U-MUSCL reconstructed variables ---*/ - const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, umusclRamp); - const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, umusclRamp); - - /// TODO: Customize the limiter function. - const Double lim_i = (delta_ij_2 + proj_i*delta_ij) / (pow(proj_i,2) + delta_ij_2); - const Double lim_j = (delta_ij_2 + proj_j*delta_ij) / (pow(proj_j,2) + delta_ij_2); - - /*--- Apply reconstruction: V_L = V_i + 0.5 * lim * dV_ij^kap ---*/ - V.i.all(iVar) += 0.5 * lim_i * proj_i; - V.j.all(iVar) -= 0.5 * lim_j * proj_j; - } -} - /*! * \brief Retrieve primitive variables for points i/j, reconstructing them if needed. * \note Density and enthalpy are recomputed from ideal gas EOS. @@ -219,17 +76,7 @@ FORCEINLINE CPair reconstructPrimitives(const Int& iEdge, if (muscl) { /*--- Reconstruct density and enthalpy without using their gradients. ---*/ constexpr auto nVarGrad = ReconVarType::nVar - 2; - switch (limiterType) { - case LIMITER::NONE: - musclUnlimited(iPoint, jPoint, vector_ij, gradients, V, kappa, umusclRamp); - break; - case LIMITER::VAN_ALBADA_EDGE: - musclEdgeLimited(iPoint, jPoint, vector_ij, gradients, V, kappa, umusclRamp); - break; - default: - musclPointLimited(iPoint, jPoint, vector_ij, limiters, gradients, V, kappa, umusclRamp); - break; - } + reconstruct(iPoint, jPoint, vector_ij, gradients, limiters, limiterType, 0, V, kappa, umusclRamp); /*--- Recompute density using the reconstructed pressure and temperature. ---*/ V.i.density() = V.i.pressure() / (gasConst * V.i.temperature()); V.j.density() = V.j.pressure() / (gasConst * V.j.temperature()); diff --git a/SU2_CFD/include/numerics_simd/util.hpp b/SU2_CFD/include/numerics_simd/util.hpp index 79594268be8..9db7721803f 100644 --- a/SU2_CFD/include/numerics_simd/util.hpp +++ b/SU2_CFD/include/numerics_simd/util.hpp @@ -1,6 +1,6 @@ /*! * \file util.hpp - * \brief Generic auxiliary functions. + * \brief Vector, matrix and index types bound to the SIMD Double/Int of CNumericsSIMD.hpp. * \author P. Gomes * \version 8.5.0 "Harrier" * @@ -28,248 +28,10 @@ #pragma once #include "CNumericsSIMD.hpp" -#include "../../../Common/include/containers/C2DContainer.hpp" -#include "../../../Common/include/linear_algebra/CSysVector.hpp" -#include "../../../Common/include/linear_algebra/CSysMatrix.hpp" +#include "../numerics/util.hpp" -/*! - * \brief Static vector and matrix types. - * \note These should be used instead of C-style arrays. - */ -template -using Vector = C2DContainer; - -template using VectorInt = Vector; -template using VectorDbl = Vector; - -template -using Matrix = C2DContainer; - -template using MatrixInt = Matrix; -template using MatrixDbl = Matrix; - -/*! - * \brief Constexpr version of max. - */ -inline constexpr size_t Max(size_t a, size_t b) { return a>b? a : b; } - -/*! - * \brief Simple pair type for i/j variables. - */ -template -struct CPair { - T i, j; -}; - -/*! - * \brief Dot product. - */ -template -FORCEINLINE Double dot(ForwardIterator iterator, const T* ptr) { - Double sum = 0.0; - for (size_t iDim = 0; iDim < nDim; ++iDim) { - sum += *(iterator++) * ptr[iDim]; - } - return sum; -} - -/*! - * \overload Dot product. - */ -template -FORCEINLINE Double dot(ForwardIterator iterator, const VectorDbl& vector) { - return dot(iterator, vector.data()); -} - -/*! - * \overload Dot product. - */ -template -FORCEINLINE Double dot(const VectorDbl& a, const VectorDbl& b) { - return dot(a.data(), b.data()); -} - -/*! - * \brief Squared norm. - */ -template -FORCEINLINE Double squaredNorm(ForwardIterator iterator) { - Double sum = 0.0; - for (size_t iDim = 0; iDim < nDim; ++iDim) { - sum += pow(*(iterator++),2); - } - return sum; -} - -/*! - * \overload Squared norm. - */ -template -FORCEINLINE Double squaredNorm(const VectorDbl& vector) { - return squaredNorm(vector.data()); -} - -/*! - * \brief Tangential projection. - */ -template -FORCEINLINE VectorDbl tangentProjection(const MatrixDbl& tensor, - const VectorDbl& unitVector) { - VectorDbl proj; - for (size_t iDim = 0; iDim < nDim; ++iDim) - proj(iDim) = dot(tensor[iDim], unitVector); - - Double normalProj = dot(proj, unitVector); - - for (size_t iDim = 0; iDim < nDim; ++iDim) - proj(iDim) -= normalProj * unitVector(iDim); +template using VectorInt = Vector; +template using VectorDbl = Vector; - return proj; -} - -/*! - * \brief Vector norm. - */ -template -FORCEINLINE Double norm(const VectorDbl& vector) { return sqrt(squaredNorm(vector)); } - -#ifndef CODI_REVERSE_TYPE -/*! - * \brief Gather a single variable from index iPoint of a 1D container. - */ -template -FORCEINLINE Double gatherVariables(Int iPoint, const Container& vars) { - return *vars.innerIter(iPoint); -} - -/*! - * \brief Gather a vector of variables (size nVar) from row iPoint of a 2D container. - */ -template -FORCEINLINE VectorDbl gatherVariables(Int iPoint, const Container& vars) { - return vars.template get >(iPoint); -} - -/*! - * \brief Gather a matrix of variables from outer index iPoint of a 3D container. - */ -template -FORCEINLINE MatrixDbl gatherVariables(Int iPoint, const Container& vars) { - return vars.template get >(iPoint); -} -#else - -namespace { - template = 0> - FORCEINLINE const su2double& get(const Container& vars, unsigned long iPoint) { return vars(iPoint); } - - /*--- When getting 1 variable from a matrix container, we assume it is the first. ---*/ - template = 0> - FORCEINLINE const su2double& get(const Container& vars, unsigned long iPoint) { return vars(iPoint,0); } -} - -template -FORCEINLINE Double gatherVariables(Int iPoint, const Container& vars) { - Double x; - for (size_t k=0; k -FORCEINLINE VectorDbl gatherVariables(Int iPoint, const Container& vars) { - VectorDbl x; - for (size_t i=0; i -FORCEINLINE MatrixDbl gatherVariables(Int iPoint, const Container& vars) { - MatrixDbl x; - for (size_t i=0; i -FORCEINLINE void stopPreacc(VectorDbl& x) { - AD::SetPreaccOut(x, nVar, Double::Size); - AD::EndPreacc(); -} - -/*! - * \brief Distance vector, from point i to point j. - */ -template -FORCEINLINE VectorDbl distanceVector(Int iPoint, Int jPoint, - const Container& coords) { - auto coord_i = gatherVariables(iPoint, coords); - auto coord_j = gatherVariables(jPoint, coords); - VectorDbl vector_ij; - for (size_t iDim = 0; iDim < nDim; ++iDim) { - vector_ij(iDim) = coord_j(iDim) - coord_i(iDim); - } - return vector_ij; -} - -/*! - * \brief Update the matrix and right-hand-side of a linear system. - */ -template -FORCEINLINE void updateLinearSystem(Int iEdge, - Int iPoint, - Int jPoint, - bool implicit, - UpdateType updateType, - Double updateMask, - const VectorDbl& flux, - const MatrixDbl& jac_i, - const MatrixDbl& jac_j, - CSysVector& vector, - SparseMatrixType& matrix) { - if (updateType == UpdateType::COLORING) { - vector.UpdateBlocks(iPoint, jPoint, flux, updateMask); - if(implicit) { - auto wasActive = AD::BeginPassive(); - matrix.SetBlocks(iEdge, iPoint, jPoint, jac_i, jac_j, updateMask); - AD::EndPassive(wasActive); - } - } - else { - vector.SetBlock(iEdge, flux, updateMask); - if(implicit) { - auto wasActive = AD::BeginPassive(); - matrix.SetBlocks(iEdge, jac_i, jac_j, updateMask); - AD::EndPassive(wasActive); - } - } -} - -/*! - * \brief Store the (scalar) mass flux of an edge, e.g. for "bounded scalar" transport equations. - * \note No-op if "target" is null. As with CEdge's Nodes/Normal, edges within a SIMD group are - * contiguous (coloring groups are multiples of the SIMD size), so this is a plain vectorized store - * starting at iEdge[0], relying on "target" being padded to a multiple of the SIMD size. - */ -FORCEINLINE void updateEdgeMassFlux(Int iEdge, - const Double& massFlux, - su2activevector* target) { - if (target) massFlux.store(&(*target)[iEdge[0]]); -} +template using MatrixInt = Matrix; +template using MatrixDbl = Matrix; diff --git a/SU2_CFD/include/solvers/CHeatSolver.hpp b/SU2_CFD/include/solvers/CHeatSolver.hpp index a0cce12c57c..1097b2cd06f 100644 --- a/SU2_CFD/include/solvers/CHeatSolver.hpp +++ b/SU2_CFD/include/solvers/CHeatSolver.hpp @@ -93,46 +93,12 @@ class CHeatSolver final : public CScalarSolver { } } - /*! - * \brief Compute the viscous flux for the scalar equation at a particular edge. - * \param[in] iEdge - Edge for which we want to compute the flux - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] numerics - Description of the numerical method. - * \param[in] config - Definition of the particular problem. - * \note Calls a generic implementation after defining a SolverSpecificNumerics object. - */ - inline void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, const CConfig* config) override { - const CVariable* flow_nodes = flow ? solver_container[FLOW_SOL]->GetNodes() : nullptr; - - const su2double const_diffusivity = config->GetThermalDiffusivity(); - const su2double pr_turb = config->GetPrandtl_Turb(); - - su2double thermal_diffusivity_i{}, thermal_diffusivity_j{}; - - /*--- Computes the thermal diffusivity to use in the viscous numerics. ---*/ - auto compute_thermal_diffusivity = [&](unsigned long iPoint, unsigned long jPoint) { - if (flow) { - thermal_diffusivity_i = flow_nodes->GetThermalConductivity(iPoint) / flow_nodes->GetSpecificHeatCp(iPoint) + - flow_nodes->GetEddyViscosity(iPoint) / pr_turb; - thermal_diffusivity_j = flow_nodes->GetThermalConductivity(jPoint) / flow_nodes->GetSpecificHeatCp(jPoint) + - flow_nodes->GetEddyViscosity(jPoint) / pr_turb; - numerics->SetDiffusionCoeff(&thermal_diffusivity_i, &thermal_diffusivity_j); - } else { - numerics->SetDiffusionCoeff(&const_diffusivity, &const_diffusivity); - } - }; - /*--- Compute residual and Jacobians. ---*/ - Viscous_Residual_impl(compute_thermal_diffusivity, iEdge, geometry, solver_container, numerics, config); - } - public: /*! * \brief Constructor of the class. */ - CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh); + CHeatSolver(CGeometry *geometry, CConfig *config, const CSolver* flow_solver, unsigned short iMesh); /*! * \brief Restart residual and compute gradients. @@ -181,13 +147,16 @@ class CHeatSolver final : public CScalarSolver { unsigned short iMesh) override; /*! - * \brief Compute the viscous residuals for the turbulent equation. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] numerics_container - Description of the numerical method. - * \param[in] config - Definition of the particular problem. - * \param[in] iMesh - Index of the mesh in multigrid computations. - * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \brief Resolve the compile-time parameters of CScalarFlux_Heat and run one of this solver's + * boundaries through the shared boundary flux pass. + * \param[in] opt - Flags of the boundary, from one of ScalarFluxOptions' named constructors. + */ + void BoundaryFlux(CGeometry* geometry, CSolver** solver_container, CConfig* config, const ScalarFluxOptions& opt, + unsigned short val_marker); + + /*! + * \brief Diffusion for solid conduction, called unconditionally unlike Upwind_Residual. A no-op + * for a fluid zone, where diffusion was already computed together with convection. */ void Viscous_Residual(CGeometry *geometry, CSolver **solver_container, @@ -273,6 +242,34 @@ class CHeatSolver final : public CScalarSolver { CConfig *config, unsigned short val_marker) override; + /*! + * \brief Impose the far-field boundary condition. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Description of the numerical method. + * \param[in] visc_numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ + void BC_Far_Field(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) override; + + /*! + * \brief Impose the fluid interface (sliding mesh) boundary condition, via the + * CScalarFlux_Heat edge kernel. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] conv_numerics - Unused, kept only for the boundary condition dispatch. + * \param[in] visc_numerics - Unused, kept only for the boundary condition dispatch. + * \param[in] config - Definition of the particular problem. + */ + void BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CNumerics *visc_numerics, CConfig *config) override; + /*! * \brief Impose the (received) conjugate heat variables. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index 500233645bf..a8e5d2aa699 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -30,11 +30,26 @@ #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../numerics/scalar/scalar_edge_flux.hpp" #include "../variables/CScalarVariable.hpp" +#include "../variables/CEulerVariable.hpp" #include "../variables/CFlowVariable.hpp" +#include "../variables/CGhostFlowVariable.hpp" +#include "../variables/CIncEulerVariable.hpp" #include "../variables/CPrimitiveIndices.hpp" #include "CSolver.hpp" +/*! + * \brief Carries a type through a value, so a runtime branch can hand a compile-time type to a + * generic lambda (its parameter deduces as CTypeTag, and the lambda recovers T as + * decltype(tag)::type). Standing in for a C++20 template lambda, which this project's + * C++17 baseline does not have. + */ +template +struct CTypeTag { + using type = T; +}; + /*! * \brief Main class for defining a scalar solver. * \tparam VariableType - Class of variable used by the solver inheriting from this template. @@ -80,6 +95,17 @@ class CScalarSolver : public CSolver { CSysVector EdgeFluxes; /*!< \brief Flux across each edge. */ CSysVector EdgeFluxesDiff; /*!< \brief Flux difference between ij and ji for non-conservative discretisation. */ + /*--- Ghost states of the marker currently being processed by a boundary, indexed by vertex + * and sized to the largest marker; same container types as the interior ones, so the flux + * kernels read a boundary through the same accessors as an interior edge. Boundary loops run + * one marker at a time, parallel over its vertices, so the buffers are written and consumed + * before the next marker reaches them (see BoundaryFluxResidual). ---*/ + unique_ptr ghostNodes; /*!< \brief Allocated by the derived solver, whose VariableType constructor it alone knows how to call. */ + unique_ptr ghostFlowNodes; /*!< \brief Sized from the flow solver, see the constructor. */ + su2activematrix ghostNormal; /*!< \brief Outward normals, sign flipped from the vertex normals. */ + su2activematrix ghostCoord; /*!< \brief Reflected coordinates, read by the diffusion sites. */ + su2vector ghostSkip; /*!< \brief Whether a vertex contributes no flux, set by the fill pass. */ + /*! * \brief The highest level in the variable hierarchy this solver can safely use. */ @@ -90,159 +116,6 @@ class CScalarSolver : public CSolver { */ inline CVariable* GetBaseClassPointerToNodes() final { return nodes; } - /*! - * \brief Compute the viscous flux for the scalar equation at a particular edge. - * \tparam SolverSpecificNumericsFunc - lambda-function, that implements solver specific contributions to numerics. - * \note The functor has to implement (iPoint, jPoint) - * \param[in] iEdge - Edge for which we want to compute the flux - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] numerics - Description of the numerical method. - * \param[in] config - Definition of the particular problem. - */ - template - FORCEINLINE void Viscous_Residual_impl(const SolverSpecificNumericsFunc& SolverSpecificNumerics, const unsigned long iEdge, - const CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, - const CConfig* config) { - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - CFlowVariable* flowNodes = solver_container[FLOW_SOL] ? - su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()) : nullptr; - - /*--- Points in edge ---*/ - - auto iPoint = geometry->edges->GetNode(iEdge, 0); - auto jPoint = geometry->edges->GetNode(iEdge, 1); - - /*--- Points coordinates, and normal vector ---*/ - - numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint)); - numerics->SetNormal(geometry->edges->GetNormal(iEdge)); - - /*--- Conservative variables w/o reconstruction ---*/ - - if (flowNodes) { - numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), flowNodes->GetPrimitive(jPoint)); - } - - /*--- Turbulent variables w/o reconstruction, and its gradients ---*/ - - numerics->SetScalarVar(nodes->GetSolution(iPoint), nodes->GetSolution(jPoint)); - numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(jPoint)); - - /*--- Call Numerics contribution which are Solver-Specifc. Implemented in the caller: Viscous_Residual. ---*/ - - SolverSpecificNumerics(iPoint, jPoint); - - /*--- Compute residual, and Jacobians ---*/ - - auto residual = numerics->ComputeResidual(config); - - if (ReducerStrategy) { - EdgeFluxes.SubtractBlock(iEdge, residual); - if (implicit) Jacobian.UpdateBlocksSub(iEdge, residual.jacobian_i, residual.jacobian_j); - } else { - LinSysRes.SubtractBlock(iPoint, residual); - LinSysRes.AddBlock(jPoint, residual); - if (implicit) Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); - } - } - - /*! - * \brief Compute the viscous flux for the turbulence equations at a particular edge for a non-conservative discretisation. - * \tparam SolverSpecificNumericsTemp - lambda-function, to implement solver specific contributions to numerics. - * \note The functor has to implement (iPoint, jPoint) - * \param[in] iEdge - Edge for which we want to compute the flux - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] numerics - Description of the numerical method. - * \param[in] config - Definition of the particular problem. - */ - template - void Viscous_Residual_NonCons(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, const CConfig* config, SolverSpecificNumericsFunc&& SolverSpecificNumerics) { - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - CFlowVariable* flowNodes = solver_container[FLOW_SOL] ? - su2staticcast_p(solver_container[FLOW_SOL]->GetNodes()) : nullptr; - - const auto iPoint = geometry->edges->GetNode(iEdge, 0); - const auto jPoint = geometry->edges->GetNode(iEdge, 1); - - /*--- Lambda function to compute the flux ---*/ - auto ComputeFlux = [&](unsigned long iPoint, unsigned long jPoint, const su2double* normal) { - numerics->SetCoord(geometry->nodes->GetCoord(iPoint),geometry->nodes->GetCoord(jPoint)); - numerics->SetNormal(normal); - - if (flowNodes) { - numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), flowNodes->GetPrimitive(jPoint)); - } - - /*--- Solver specific numerics contribution. ---*/ - SolverSpecificNumerics(iPoint, jPoint); - - numerics->SetScalarVar(nodes->GetSolution(iPoint), nodes->GetSolution(jPoint)); - numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(jPoint)); - - return numerics->ComputeResidual(config); - }; - - /*--- Compute fluxes and jacobians i->j ---*/ - const su2double* normal = geometry->edges->GetNormal(iEdge); - auto residual_ij = ComputeFlux(iPoint, jPoint, normal); - - su2mixedfloat *Block_ii = nullptr, *Block_ij = nullptr, *Block_ji = nullptr, *Block_jj = nullptr; - if (implicit) { - Jacobian.GetBlocks(iEdge, iPoint, jPoint, Block_ii, Block_ij, Block_ji, Block_jj); - } - if (ReducerStrategy) { - EdgeFluxes.SubtractBlock(iEdge, residual_ij); - EdgeFluxesDiff.SetBlock(iEdge, residual_ij); - if (implicit) { - /*--- For the reducer strategy the Jacobians are averaged for simplicity. ---*/ - for (int iVar=0; iVari ---*/ - su2double flipped_normal[MAXNDIM]; - for (auto iDim = 0u; iDim < nDim; iDim++) flipped_normal[iDim] = -normal[iDim]; - - auto residual_ji = ComputeFlux(jPoint, iPoint, flipped_normal); - if (ReducerStrategy) { - EdgeFluxesDiff.AddBlock(iEdge, residual_ji); - if (implicit) { - for (int iVar=0; iVar. + * \param[in] opt - Loop invariant flags built by the caller from the current CConfig state. */ - inline virtual void Viscous_Residual(const unsigned long iEdge, const CGeometry* geometry, CSolver** solver_container, - CNumerics* numerics, const CConfig* config) { - /*--- Define an empty object for solver specific numerics contribution. In case there are none, this default - *--- implementation will be called ---*/ - auto SolverSpecificNumerics = [&](unsigned long iPoint, unsigned long jPoint) {}; + template + void EdgeFluxResidual(const CGeometry* geometry, CSolver** solver_container, const CConfig* config, + const ScalarFluxOptions& opt); - /*--- Now instantiate the generic implementation with the functor above. ---*/ + /*! + * \brief Write the four flow primitives the flux kernels read into one row of ghostFlowNodes. + * \param[in] iVertex - Vertex of the marker currently being processed. + * \param[in] V - Row of flow primitives to copy from (e.g. GetCharacPrimVar's or a sliding state's). + */ + inline void SetGhostPrimitives(unsigned long iVertex, const su2double* V) { + auto* ghostV = ghostFlowNodes->GetPrimitive(iVertex); + ghostV[prim_idx.Density()] = V[prim_idx.Density()]; + for (auto iDim = 0u; iDim < nDim; ++iDim) ghostV[prim_idx.Velocity() + iDim] = V[prim_idx.Velocity() + iDim]; + ghostV[prim_idx.LaminarViscosity()] = V[prim_idx.LaminarViscosity()]; + ghostV[prim_idx.EddyViscosity()] = V[prim_idx.EddyViscosity()]; + /*--- NEMO's primitive layout has no single thermal conductivity or specific heat, so its + * CIndices returns a sentinel for these two; only heat reads them, and NEMO rejects any + * scalar transport at configuration. ---*/ + if (prim_idx.ThermalConductivity() != std::numeric_limits::max()) { + ghostV[prim_idx.ThermalConductivity()] = V[prim_idx.ThermalConductivity()]; + ghostV[prim_idx.CpTotal()] = V[prim_idx.CpTotal()]; + } + } - Viscous_Residual_impl(SolverSpecificNumerics, iEdge, geometry, solver_container, numerics, config); + /*! + * \brief Generic boundary flux pass, run after a boundary's fill pass has written the ghost + * row, the outward normal and (for the diffusion sites) the ghost gradient of every + * vertex of the marker. The ghost point has no row, so only the contribution to the + * interior point is assembled. + * \tparam Scheme - Same model the interior loop uses, instantiated with muscl false. + */ + template + void BoundaryFluxResidual(const CGeometry* geometry, CSolver** solver_container, const CConfig* config, + const ScalarFluxOptions& opt, unsigned short val_marker); + + /*! + * \brief Generic fluid interface (sliding mesh) flux pass, shared by every model. The convective + * term is a per-donor weighted average, computed in the same pass that fills the ghost row + * of each donor; the diffusive term is computed once per vertex, after the donor loop, + * from the ghost state the last donor left behind. This does not fit the + * fill-pass-then-BoundaryFluxResidual shape the other boundaries use, so it drives the + * kernel directly. + * \tparam Scheme - Same model the interior loop uses, instantiated with muscl false. + * \param[in] fillGhostExtras - Functor (iVertex, iPoint) writing the auxiliary ghost fields the + * model's diffusion coefficients read, e.g. SST's blending function or the species + * mass diffusivities. Called once per vertex, before the diffusive flux. + */ + template + void FluidInterfaceFluxResidual(const CGeometry* geometry, CSolver** solver_container, const CConfig* config, + const ScalarFluxOptions& optConv, const ScalarFluxOptions& optVisc, + const GhostFunc& fillGhostExtras); + + /*! + * \brief Write the outward normal of one vertex into the ghost row and mark the vertex as + * contributing a flux. + * \note Vertex normals point into the domain, the flux convention needs them outward. + */ + inline void SetGhostGeometry(const CGeometry* geometry, unsigned short val_marker, unsigned long iVertex) { + for (auto iDim = 0u; iDim < nDim; ++iDim) + ghostNormal(iVertex, iDim) = -geometry->vertex[val_marker][iVertex]->GetNormal(iDim); + ghostSkip[iVertex] = false; + } + + /*! + * \brief Write what the diffusion term of a boundary reads beyond the ghost solution: the + * coordinate of the interior point reflected about the boundary, and the interior + * gradient mirrored into the ghost row. + * \param[in] iPoint - Interior point of the vertex. + * \param[in] jPoint - Point the interior one is reflected about, the vertex's normal neighbor. + */ + inline void SetGhostDiffusionState(const CGeometry* geometry, unsigned long iVertex, unsigned long iPoint, + unsigned long jPoint) { + su2double Coord_Reflected[MAXNDIM]; + GeometryToolbox::PointPointReflect(nDim, geometry->nodes->GetCoord(jPoint), geometry->nodes->GetCoord(iPoint), + Coord_Reflected); + for (auto iDim = 0u; iDim < nDim; ++iDim) ghostCoord(iVertex, iDim) = Coord_Reflected[iDim]; + + auto ghostGrad = ghostNodes->GetGradient(iVertex); + const auto interiorGrad = nodes->GetGradient(iPoint); + for (auto iVar = 0u; iVar < nVar; ++iVar) + for (auto iDim = 0u; iDim < nDim; ++iDim) ghostGrad(iVar, iDim) = interiorGrad(iVar, iDim); + } + + /*! + * \brief Resolve the compile-time parameters of a scalar flux kernel, the flow indices, the + * dimension and the equation count, from the runtime state, and call f with a CTypeTag of + * the resulting scheme type: f is a generic lambda, + * `[&](auto tag){ using Scheme = typename decltype(tag)::type; ... }`. + * \tparam Model - Model class template, e.g. CScalarFlux_SST, taking the four parameters of + * CUpwScalarBase: value type, flow indices, dimension and equation count. + * \tparam nVarList - Equation counts to instantiate. Dynamic matches any count, a static one is + * taken when it equals the solver's nVar; the counts are tried in the order given. + * \note NEMO is not one of the index branches: transported scalars are rejected for it at + * configuration. + */ + template