Skip to content

Fix streamwise periodic heat, massflow and auxiliary gradients #2965

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
rois1995 wants to merge 11 commits into
base: develop
Choose a base branch
from
Draft

Fix streamwise periodic heat, massflow and auxiliary gradients #2965

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
e4deaa2
Stop when streamwise periodicity is used with multigrid
Oct 5, 2026
071ad4d
Use INTEGRATED_HEATFLUX in the integrated heat of streamwise periodic…
Oct 5, 2026
8982549
Fix the turbulent source term of the streamwise periodic temperature
Oct 5, 2026
a7a424b
Use the signed distance to recover the streamwise periodic pressure a…
Oct 5, 2026
192b327
Initialize the streamwise periodic values
Oct 5, 2026
0c8e453
Update streamwise periodic pressure drop once per OpenMP team
Oct 7, 2026
609d207
Complete periodic Green-Gauss auxiliary gradients
Oct 7, 2026
cbf1003
Exchange auxiliary periodic gradients and correct streamwise heated-w…
Oct 7, 2026
ece5279
Normalize assignment spacing reported by CodeFactor
Oct 7, 2026
95c84b7
Promote periodic auxiliary matrix count before multiplication
Oct 7, 2026
4d62254
Merge develop into fix_streamwise_periodic
Oct 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 8 additions & 6 deletions Common/include/option_structure.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2865,6 +2865,8 @@ enum PERIODIC_QUANTITIES {
PERIODIC_LIM_PRIM_1 , /*!< \brief Primitive limiter communication phase 1 of 2 (periodic only). */
PERIODIC_LIM_PRIM_2 , /*!< \brief Primitive limiter communication phase 2 of 2 (periodic only). */
PERIODIC_IMPLICIT , /*!< \brief Implicit update communication to ensure consistency across periodic boundaries. */
PERIODIC_AUXVAR_LS , /*!< \brief Auxiliary scalar gradient communication for weighted least-squares. */
PERIODIC_AUXVAR_GG , /*!< \brief Auxiliary scalar gradient communication for Green-Gauss. */
};

/*!
Expand Down Expand Up @@ -3001,12 +3003,12 @@ static const MapType<std::string, ENUM_STREAMWISE_PERIODIC> Streamwise_Periodic_
* \brief Container to hold Variables for streamwise Periodic flow as they are often used together in places.
*/
struct StreamwisePeriodicValues {
su2double Streamwise_Periodic_PressureDrop; /*!< \brief Value of prescribed pressure drop [Pa] which results in an artificial body force vector. */
su2double Streamwise_Periodic_MassFlow; /*!< \brief Value of current massflow [kg/s] which results in a delta p and therefore an artificial body force vector. */
su2double Streamwise_Periodic_IntegratedHeatFlow; /*!< \brief Value of of the net sum of heatflow [W] into the domain. */
su2double Streamwise_Periodic_InletTemperature; /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */
su2double Streamwise_Periodic_BoundaryArea; /*!< \brief Global Surface area of the streamwise periodic interface. */
su2double Streamwise_Periodic_AvgDensity; /*!< \brief Area avg density on the periodic interface. */
su2double Streamwise_Periodic_PressureDrop = 0.0; /*!< \brief Value of prescribed pressure drop [Pa] which results in an artificial body force vector. */
su2double Streamwise_Periodic_MassFlow = 0.0; /*!< \brief Value of current massflow [kg/s] which results in a delta p and therefore an artificial body force vector. */
su2double Streamwise_Periodic_IntegratedHeatFlow = 0.0; /*!< \brief Value of of the net sum of heatflow [W] into the domain. */
su2double Streamwise_Periodic_InletTemperature = 0.0; /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */
su2double Streamwise_Periodic_BoundaryArea = 0.0; /*!< \brief Global Surface area of the streamwise periodic interface. */
su2double Streamwise_Periodic_AvgDensity = 0.0; /*!< \brief Area avg density on the periodic interface. */
};

/*!
Expand Down
17 changes: 17 additions & 0 deletions Common/src/CConfig.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@

#include "../include/basic_types/ad_structure.hpp"
#include "../include/toolboxes/printing_toolbox.hpp"
#include "../include/toolboxes/geometry_toolbox.hpp"
#include "../include/toolboxes/SwapBytes.hpp"

using namespace PrintingToolbox;
Expand Down Expand Up @@ -5558,12 +5559,28 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
SU2_MPI::Error("Streamwise Periodic Flow + Incompressible Euler: Not tested yet.", CURRENT_FUNCTION);
if (nMarker_PerBound == 0)
SU2_MPI::Error("A MARKER_PERIODIC pair has to be set with KIND_STREAMWISE_PERIODIC != NONE.", CURRENT_FUNCTION);
/*--- The first pair defines the streamwise direction used by every source. ---*/
const auto* translation = GetPeriodic_Translation(0);
if (GeometryToolbox::SquaredNorm(3, translation) == 0)
SU2_MPI::Error("The first streamwise periodic pair must have a nonzero translation.", CURRENT_FUNCTION);
const auto* angles = GetPeriodicRotAngles(Marker_PerBound[0]);
su2double rotation[3][3];
GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation);
for (auto i = 0u; i < 3; ++i)
for (auto j = 0u; j < 3; ++j)
if (fabs(rotation[i][j] - (i == j)) > 1e-12)
SU2_MPI::Error("The first streamwise periodic pair must be purely translational.", CURRENT_FUNCTION);
if (Energy_Equation && Streamwise_Periodic_Temperature && (nMarker_HeatTransfer || nMarker_CHTInterface))
SU2_MPI::Error("STREAMWISE_PERIODIC_TEMPERATURE requires prescribed MARKER_HEATFLUX walls; "
"MARKER_HEATTRANSFER and MARKER_CHT_INTERFACE are unsupported.", CURRENT_FUNCTION);
if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal != 0)
SU2_MPI::Error("No MARKER_ISOTHERMAL marker allowed with STREAMWISE_PERIODIC_TEMPERATURE= YES, only MARKER_HEATFLUX & MARKER_SYM.", CURRENT_FUNCTION);
if (Ref_Inc_NonDim != DIMENSIONAL)
SU2_MPI::Error("Streamwise Periodicity only works with \"INC_NONDIM= DIMENSIONAL\", the nondimensionalization with source terms doesn;t work in general.", CURRENT_FUNCTION);
if (Axisymmetric)
SU2_MPI::Error("Streamwise Periodicity terms does not not have axisymmetric corrections.", CURRENT_FUNCTION);
if (nMGLevels != 0)
SU2_MPI::Error("Streamwise Periodicity does not support multigrid, use MGLEVEL= 0.", CURRENT_FUNCTION);
if (!Energy_Equation) Streamwise_Periodic_Temperature = false;
} else {
/*--- Safety measure ---*/
Expand Down
6 changes: 4 additions & 2 deletions SU2_CFD/src/numerics/flow/flow_sources.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -794,8 +794,10 @@ CNumerics::ResidualType<> CSourceIncStreamwise_Periodic::ComputeResidual(const C
/*--- If a RANS turbulence model is used, an additional source term, based on the eddy viscosity gradient is added. ---*/
if(turbulent) {

/*--- Compute a scalar factor ---*/
scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * sqrt(norm2_translation) * Prandtl_Turb);
/*--- Compute a scalar factor. The temperature is T = T_periodic + Q / (mdot * cp * |t|^2) * dot_prod(t*x)
(Patankar et al., J. Heat Transfer 99(2), 1977), so the turbulent heat flux cp * mu_t / Pr_t * grad(T)
has the extra divergence Q / (mdot * |t|^2 * Pr_t) * dot_prod(t*grad(mu_t)). ---*/
scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * norm2_translation * Prandtl_Turb);

/*--- Compute scalar product between periodic translation vector and eddy viscosity gradient. ---*/
dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, AuxVar_Grad_i[0]);
Expand Down
30 changes: 17 additions & 13 deletions SU2_CFD/src/solvers/CIncEulerSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1889,19 +1889,23 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
/*--- The Pressure drop is iteratively adapted to result in the prescribed Target-Massflow. ---*/
/*---------------------------------------------------------------------------------------------*/

/*--- Compute update to Delta p based on massflow-difference ---*/
const su2double Average_Density_Global = SPvals.Streamwise_Periodic_AvgDensity;
const su2double Area_Global = SPvals.Streamwise_Periodic_BoundaryArea;
const su2double TargetMassFlow = config->GetStreamwise_Periodic_TargetMassFlow() / (config->GetDensity_Ref() * config->GetVelocity_Ref());
const su2double MassFlow_Global = SPvals.Streamwise_Periodic_MassFlow;
const su2double ddP = 0.5 / ( Average_Density_Global * pow(Area_Global, 2)) * (pow(TargetMassFlow, 2) - pow(MassFlow_Global, 2));

/*--- Store updated pressure difference ---*/
const su2double damping_factor = config->GetInc_Outlet_Damping();
SPvalsUpdated = SPvals;
SPvalsUpdated.Streamwise_Periodic_PressureDrop += damping_factor*ddP;
if (!config->GetDiscrete_Adjoint())
SPvals = SPvalsUpdated;
BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
{
/*--- Compute update to Delta p based on massflow-difference ---*/
const su2double Average_Density_Global = SPvals.Streamwise_Periodic_AvgDensity;
const su2double Area_Global = SPvals.Streamwise_Periodic_BoundaryArea;
const su2double TargetMassFlow = config->GetStreamwise_Periodic_TargetMassFlow() / (config->GetDensity_Ref() * config->GetVelocity_Ref());
const su2double MassFlow_Global = SPvals.Streamwise_Periodic_MassFlow;
const su2double ddP = 0.5 / ( Average_Density_Global * pow(Area_Global, 2)) * (pow(TargetMassFlow, 2) - pow(MassFlow_Global, 2));

/*--- Store updated pressure difference ---*/
const su2double damping_factor = config->GetInc_Outlet_Damping();
SPvalsUpdated = SPvals;
SPvalsUpdated.Streamwise_Periodic_PressureDrop += damping_factor*ddP;
if (!config->GetDiscrete_Adjoint())
SPvals = SPvalsUpdated;
}
END_SU2_OMP_SAFE_GLOBAL_ACCESS

/*--- Set delta_p, m_dot, inlet_T, integrated_heat ---*/
numerics->SetStreamwisePeriodicValues(SPvalsUpdated);
Expand Down
10 changes: 6 additions & 4 deletions SU2_CFD/src/solvers/CIncNSSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,8 @@ void CIncNSSolver::GetStreamwise_Periodic_Properties(const CGeometry *geometry,

/*--- Identify the boundary by string name and retrive heatflux from config ---*/
const auto Marker_StringTag = config->GetMarker_All_TagBound(iMarker);
const su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_StringTag);
su2double Wall_HeatFlux = config->GetWall_HeatFlux(Marker_StringTag);
if (config->GetIntegrated_HeatFlux()) Wall_HeatFlux /= geometry->GetSurfaceArea(config, iMarker);

for (auto iVertex = 0ul; iVertex < geometry->nVertex[iMarker]; iVertex++) {

Expand Down Expand Up @@ -254,8 +255,8 @@ void CIncNSSolver::Compute_Streamwise_Periodic_Recovered_Values(CConfig *config,

/*--- First, compute helping terms based on relative distance (0,l) between periodic markers ---*/
su2double dot_product = 0.0;
for (unsigned short iDim = 0; iDim < nDim; iDim++)
dot_product += fabs( (geometry->nodes->GetCoord(iPoint,iDim) - ReferenceNode[iDim]) * config->GetPeriodic_Translation(0)[iDim]);
for (auto iDim = 0u; iDim < nDim; iDim++)
dot_product += (geometry->nodes->GetCoord(iPoint,iDim) - ReferenceNode[iDim]) * config->GetPeriodic_Translation(0)[iDim];

/*--- Second, substract/add correction from reduced pressure/temperature to get recoverd pressure/temperature ---*/
const su2double Pressure_Recovered = nodes->GetPressure(iPoint) - SPvals.Streamwise_Periodic_PressureDrop /
Expand Down Expand Up @@ -537,7 +538,8 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
/*--- Dot product ---*/
const su2double dot_product = GeometryToolbox::DotProduct(nDim, config->GetPeriodic_Translation(0), Normal);

LinSysRes(iPoint, nDim+1) += scalar_factor*dot_product;
/*--- Normal points into the fluid: q_periodic = q + k grad(T_recovered).n. ---*/
LinSysRes(iPoint, nDim+1) -= scalar_factor*dot_product;
} // if streamwise_periodic
break;

Expand Down
42 changes: 33 additions & 9 deletions SU2_CFD/src/solvers/CSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -229,6 +229,12 @@
COUNT_PER_POINT = 2;
MPI_TYPE = COMM_TYPE::DOUBLE;
break;
case PERIODIC_AUXVAR_GG:
COUNT_PER_POINT = base_nodes->GetnAuxVar()*nDim;
MPI_TYPE = COMM_TYPE::DOUBLE;
ICOUNT = base_nodes->GetnAuxVar();
JCOUNT = nDim;
break;
case PERIODIC_SOL_GG:
case PERIODIC_SOL_GG_R:
COUNT_PER_POINT = nVar*nDim;
Expand All @@ -243,6 +249,12 @@
ICOUNT = nPrimVarGrad;
JCOUNT = nDim;
break;
case PERIODIC_AUXVAR_LS:
COUNT_PER_POINT = static_cast<unsigned long>(nDim)*nDim + base_nodes->GetnAuxVar()*nDim;
MPI_TYPE = COMM_TYPE::DOUBLE;
ICOUNT = base_nodes->GetnAuxVar();
JCOUNT = nDim;
break;
case PERIODIC_SOL_LS:
case PERIODIC_SOL_ULS:
case PERIODIC_SOL_LS_R:
Expand Down Expand Up @@ -296,6 +308,9 @@
case PERIODIC_PRIM_ULS:
return nodes->GetGradient_Primitive();
break;
case PERIODIC_AUXVAR_LS:
case PERIODIC_AUXVAR_GG:
return nodes->GetAuxVarGradient();
case PERIODIC_SOL_GG:
case PERIODIC_SOL_LS:
case PERIODIC_SOL_ULS:
Expand All @@ -309,6 +324,9 @@

const su2activematrix& selectField(CVariable* nodes, unsigned short commType) {
switch(commType) {
case PERIODIC_AUXVAR_LS:
case PERIODIC_AUXVAR_GG:
return nodes->GetAuxVar();
case PERIODIC_PRIM_GG:
case PERIODIC_PRIM_LS:
case PERIODIC_PRIM_ULS:
Expand Down Expand Up @@ -356,6 +374,8 @@

bool boundary_i, boundary_j;
bool weighted = true;
const bool auxiliary = commType == PERIODIC_AUXVAR_LS;
const auto nRotVars = max<unsigned long>(max(nVar, nPrimVar), base_nodes->GetnAuxVar());

unsigned short iVar, jVar, iDim;
unsigned short nNeighbor = 0;
Expand All @@ -372,8 +392,8 @@
auto *Und_Lapl = new su2double[nVar];
auto *Sol_Min = new su2double[std::max(nVar, nPrimVarGrad)];
auto *Sol_Max = new su2double[std::max(nVar, nPrimVarGrad)];
auto *rotPrim_i = new su2double[std::max(nVar, nPrimVar)];
auto *rotPrim_j = new su2double[std::max(nVar, nPrimVar)];
auto *rotPrim_i = new su2double[nRotVars];
auto *rotPrim_j = new su2double[nRotVars];

su2double Sensor_i = 0.0, Sensor_j = 0.0, Pressure_i, Pressure_j;
const su2double *Coord_i, *Coord_j;
Expand Down Expand Up @@ -486,499 +506,501 @@
/*--- Load the send buffers depending on the particular value
that has been requested for communication. ---*/

switch (commType) {

case PERIODIC_VOLUME:

/*--- Load the volume of the current periodic CV so that
we can accumulate the total control volume size on all
periodic faces. ---*/

bufDSend[buf_offset] = geometry->nodes->GetVolume(iPoint) +
geometry->nodes->GetPeriodicVolume(iPoint);

break;

case PERIODIC_NEIGHBORS:

nNeighbor = 0;
for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {

/*--- Check if this neighbor lies on the periodic face so
that we avoid double counting neighbors on both sides. If
not, increment the count of neighbors for the donor. ---*/

if (!geometry->nodes->GetPeriodicBoundary(jPoint))
nNeighbor++;
}

/*--- Store the number of neighbors in bufffer. ---*/

bufSSend[buf_offset] = nNeighbor;

break;

case PERIODIC_RESIDUAL:

/*--- Communicate the residual from our partial control
volume to the other side of the periodic face. ---*/

for (iVar = 0; iVar < nVar; iVar++) {
bufDSend[buf_offset+iVar] = LinSysRes(iPoint, iVar);
}

/*--- Rotate the momentum components of the residual array. ---*/

if (rotate_periodic) {
Rotate(zeros, &LinSysRes(iPoint,1), &bufDSend[buf_offset+1]);
}
buf_offset += nVar;

/*--- Load the time step for the current point. ---*/

bufDSend[buf_offset] = base_nodes->GetDelta_Time(iPoint);
buf_offset++;

/*--- For implicit calculations, we will communicate the
contributions to the Jacobian block diagonal, i.e., the
impact of the point upon itself, J_ii. ---*/

if (implicit_periodic) {

const auto block = Jacobian.GetBlockView(iPoint, iPoint);

for (iVar = 0; iVar < nVar; iVar++) {
for (jVar = 0; jVar < nVar; jVar++) {
jacBlock[iVar][jVar] = block(iVar, jVar);
}
}

/*--- Rotate the momentum columns of the Jacobian. ---*/

if (rotate_periodic) {
for (iVar = 0; iVar < nVar; iVar++) {
if (nDim == 2) {
jacBlock[1][iVar] = rotMatrix2D[0][0]*block(1, iVar) + rotMatrix2D[0][1]*block(2, iVar);
jacBlock[2][iVar] = rotMatrix2D[1][0]*block(1, iVar) + rotMatrix2D[1][1]*block(2, iVar);
} else {
jacBlock[1][iVar] = rotMatrix3D[0][0]*block(1, iVar) + rotMatrix3D[0][1]*block(2, iVar) +
rotMatrix3D[0][2]*block(3, iVar);
jacBlock[2][iVar] = rotMatrix3D[1][0]*block(1, iVar) + rotMatrix3D[1][1]*block(2, iVar) +
rotMatrix3D[1][2]*block(3, iVar);
jacBlock[3][iVar] = rotMatrix3D[2][0]*block(1, iVar) + rotMatrix3D[2][1]*block(2, iVar) +
rotMatrix3D[2][2]*block(3, iVar);
}
}
}

/*--- Load the Jacobian terms into the buffer for sending. ---*/

for (iVar = 0; iVar < nVar; iVar++) {
for (jVar = 0; jVar < nVar; jVar++) {
bufDSend[buf_offset] = jacBlock[iVar][jVar];
buf_offset++;
}
}
}

break;

case PERIODIC_IMPLICIT:

/*--- Communicate the solution from our master set of periodic
nodes (from the linear solver perspective) to the passive
periodic nodes on the matching face. This is done at the
end of the iteration to synchronize the solution after the
linear solve. ---*/

for (iVar = 0; iVar < nVar; iVar++) {
bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar);
}

/*--- Rotate the momentum components of the solution array. ---*/

if (rotate_periodic) {
Rotate(zeros, &base_nodes->GetSolution(iPoint)[1], &bufDSend[buf_offset+1]);
}

break;

case PERIODIC_LAPLACIAN:

/*--- For JST, the undivided Laplacian must be computed
consistently by using the complete control volume info
from both sides of the periodic face. ---*/

for (iVar = 0; iVar < nVar; iVar++)
Und_Lapl[iVar] = 0.0;

for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {

/*--- Avoid periodic boundary points so that we do not
duplicate edges on both sides of the periodic BC. ---*/

if (!geometry->nodes->GetPeriodicBoundary(jPoint)) {

/*--- Solution differences ---*/

for (iVar = 0; iVar < nVar; iVar++)
Diff[iVar] = (base_nodes->GetSolution(iPoint, iVar) -
base_nodes->GetSolution(jPoint,iVar));

/*--- Correction for compressible flows (use enthalpy) ---*/

if (!(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE)) {
Pressure_i = base_nodes->GetPressure(iPoint);
Pressure_j = base_nodes->GetPressure(jPoint);
Diff[nVar-1] = ((base_nodes->GetSolution(iPoint,nVar-1) + Pressure_i) -
(base_nodes->GetSolution(jPoint,nVar-1) + Pressure_j));
}

boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint);
boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint);

/*--- Both points inside the domain, or both in the boundary ---*/
/*--- iPoint inside the domain, jPoint on the boundary ---*/

if (!boundary_i || boundary_j) {
if (geometry->nodes->GetDomain(iPoint)){
for (iVar = 0; iVar< nVar; iVar++)
Und_Lapl[iVar] -= Diff[iVar];
}
}
}
}

/*--- Store the components to be communicated in the buffer. ---*/

for (iVar = 0; iVar < nVar; iVar++)
bufDSend[buf_offset+iVar] = Und_Lapl[iVar];

/*--- Rotate the momentum components of the Laplacian. ---*/

if (rotate_periodic) {
Rotate(zeros, &Und_Lapl[1], &bufDSend[buf_offset+1]);
}

break;

case PERIODIC_MAX_EIG:

/*--- Simple summation of eig calc on both periodic faces. ---*/

bufDSend[buf_offset] = base_nodes->GetLambda(iPoint);

break;

case PERIODIC_SENSOR: {
const bool msw = config->GetKind_Upwind_Flow() == UPWIND::MSW;

/*--- For the centered schemes, the sensor must be computed
consistently using info from the entire control volume
on both sides of the periodic face. ---*/

Sensor_i = 0; Sensor_j = 0;
for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {

/*--- Avoid halos and boundary points so that we don't
duplicate edges on both sides of the periodic BC. ---*/

if (geometry->nodes->GetPeriodicBoundary(jPoint)) continue;

/*--- Use density instead of pressure for incomp. flows. ---*/

if (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) {
Pressure_i = base_nodes->GetDensity(iPoint);
Pressure_j = base_nodes->GetDensity(jPoint);
} else {
Pressure_i = base_nodes->GetPressure(iPoint);
Pressure_j = base_nodes->GetPressure(jPoint);
}

boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint);
boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint);

/*--- Both points inside domain, or both on boundary ---*/
/*--- iPoint inside the domain, jPoint on the boundary ---*/

if ((!boundary_i || boundary_j) && geometry->nodes->GetDomain(iPoint)) {
if (msw) {
Sensor_i = fmax(Sensor_i, fabs(Pressure_j - Pressure_i)) / fmin(Pressure_i, Pressure_j);
} else {
Sensor_i += (Pressure_j - Pressure_i);
Sensor_j += (Pressure_i + Pressure_j);
}
}

}

/*--- Store the sensor increments to buffer. After summing
all contributions, these will be divided. ---*/

bufDSend[buf_offset] = Sensor_i;
buf_offset++;
bufDSend[buf_offset] = Sensor_j;

} break;

case PERIODIC_SOL_GG:
case PERIODIC_SOL_GG_R:
case PERIODIC_PRIM_GG:
case PERIODIC_PRIM_GG_R:
case PERIODIC_AUXVAR_GG:

/*--- Access and rotate the partial G-G gradient. These will be
summed on both sides of the periodic faces before dividing
by the volume to complete the Green-Gauss gradient calc. ---*/

for (iVar = 0; iVar < ICOUNT; iVar++) {
for (iDim = 0; iDim < nDim; iDim++) {
jacBlock[iVar][iDim] = gradient(iPoint, iVar, iDim);
}
}

/*--- Rotate the gradients in x,y,z space for all variables. ---*/

for (iVar = 0; iVar < ICOUNT; iVar++) {
Rotate(zeros, jacBlock[iVar], rotBlock[iVar]);
}

/*--- Rotate the vector components of the solution. ---*/
/*--- Auxiliary fields are scalars: rotate their spatial gradients only. ---*/

if (rotate_periodic) {
if (rotate_periodic && commType != PERIODIC_AUXVAR_GG) {
for (iDim = 0; iDim < nDim; iDim++) {
su2double d_diDim[3] = {0.0};
for (iVar = 1; iVar < 1+nDim; ++iVar) {
d_diDim[iVar-1] = rotBlock(iVar, iDim);
}
su2double rotated[3] = {0.0};
Rotate(zeros, d_diDim, rotated);
for (iVar = 1; iVar < 1+nDim; ++iVar) {
rotBlock(iVar, iDim) = rotated[iVar-1];
}
}
}

/*--- Store the partial gradient in the buffer. ---*/

for (iVar = 0; iVar < ICOUNT; iVar++) {
for (iDim = 0; iDim < nDim; iDim++) {
bufDSend[buf_offset+iVar*nDim+iDim] = rotBlock[iVar][iDim];
}
}

break;

case PERIODIC_AUXVAR_LS:
case PERIODIC_SOL_LS: case PERIODIC_SOL_ULS:
case PERIODIC_SOL_LS_R: case PERIODIC_SOL_ULS_R:
case PERIODIC_PRIM_LS: case PERIODIC_PRIM_ULS:
case PERIODIC_PRIM_LS_R: case PERIODIC_PRIM_ULS_R:

/*--- For L-S gradient calculations with rotational periodicity,
we will need to rotate the x,y,z components. To make the process
easier, we choose to rotate the initial periodic point and their
neighbor points into their location on the donor marker before
computing the terms that we need to communicate. ---*/

/*--- Set a flag for unweighted or weighted least-squares. ---*/

switch(commType) {
case PERIODIC_SOL_ULS:
case PERIODIC_SOL_ULS_R:
case PERIODIC_PRIM_ULS:
case PERIODIC_PRIM_ULS_R:
weighted = false;
break;
default:
weighted = true;
break;
}

/*--- Get coordinates for the current point. ---*/

Coord_i = geometry->nodes->GetCoord(iPoint);

/*--- Get the position vector from rotation center to point. ---*/

GeometryToolbox::Distance(nDim, Coord_i, center, distance);

/*--- Compute transformed point coordinates. ---*/

Rotate(translation, distance, rotCoord_i);

/*--- Get conservative solution and rotate if necessary. ---*/

for (iVar = 0; iVar < ICOUNT; iVar++)
rotPrim_i[iVar] = field(iPoint, iVar);

if (rotate_periodic) {
if (rotate_periodic && !auxiliary) {
Rotate(zeros, &field(iPoint,1), &rotPrim_i[1]);
}

/*--- Inizialization of variables ---*/

Cvector.resize(ICOUNT,nDim) = su2double(0.0);

r11 = 0.0; r12 = 0.0; r22 = 0.0;
r13 = 0.0; r23_a = 0.0; r23_b = 0.0; r33 = 0.0;

for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {

/*--- Avoid periodic boundary points so that we do not
duplicate edges on both sides of the periodic BC. ---*/

if (!geometry->nodes->GetPeriodicBoundary(jPoint)) {

/*--- Get coordinates for the neighbor point. ---*/

Coord_j = geometry->nodes->GetCoord(jPoint);

/*--- Get the position vector from rotation center. ---*/

GeometryToolbox::Distance(nDim, Coord_j, center, distance);

/*--- Compute transformed point coordinates. ---*/

Rotate(translation, distance, rotCoord_j);

/*--- Get conservative solution and rotate if necessary. ---*/

for (iVar = 0; iVar < ICOUNT; iVar++)
rotPrim_j[iVar] = field(jPoint,iVar);

if (rotate_periodic) {
if (rotate_periodic && !auxiliary) {
Rotate(zeros, &field(jPoint,1), &rotPrim_j[1]);
}

if (weighted) {
weight = GeometryToolbox::SquaredDistance(nDim, rotCoord_j, rotCoord_i);
} else {
weight = 1.0;
}

/*--- Sumations for entries of upper triangular matrix R ---*/

if (weight != 0.0) {

r11 += ((rotCoord_j[0]-rotCoord_i[0])*
(rotCoord_j[0]-rotCoord_i[0])/weight);
r12 += ((rotCoord_j[0]-rotCoord_i[0])*
(rotCoord_j[1]-rotCoord_i[1])/weight);
r22 += ((rotCoord_j[1]-rotCoord_i[1])*
(rotCoord_j[1]-rotCoord_i[1])/weight);

if (nDim == 3) {
r13 += ((rotCoord_j[0]-rotCoord_i[0])*
(rotCoord_j[2]-rotCoord_i[2])/weight);
r23_a += ((rotCoord_j[1]-rotCoord_i[1])*
(rotCoord_j[2]-rotCoord_i[2])/weight);
r23_b += ((rotCoord_j[0]-rotCoord_i[0])*
(rotCoord_j[2]-rotCoord_i[2])/weight);
r33 += ((rotCoord_j[2]-rotCoord_i[2])*
(rotCoord_j[2]-rotCoord_i[2])/weight);
}

/*--- Entries of c:= transpose(A)*b ---*/

for (iVar = 0; iVar < ICOUNT; iVar++)
for (iDim = 0; iDim < nDim; iDim++)
Cvector(iVar,iDim) += ((rotCoord_j[iDim]-rotCoord_i[iDim])*
(rotPrim_j[iVar]-rotPrim_i[iVar])/weight);

}
}
}

/*--- We store and communicate the increments for the matching
upper triangular matrix (weights) and the r.h.s. vector.
These will be accumulated before completing the L-S gradient
calculation for each periodic point. ---*/

if (nDim == 2) {
bufDSend[buf_offset] = r11; buf_offset++;
bufDSend[buf_offset] = r12; buf_offset++;
bufDSend[buf_offset] = 0.0; buf_offset++;
bufDSend[buf_offset] = r22; buf_offset++;
}
if (nDim == 3) {
bufDSend[buf_offset] = r11; buf_offset++;
bufDSend[buf_offset] = r12; buf_offset++;
bufDSend[buf_offset] = r13; buf_offset++;

bufDSend[buf_offset] = 0.0; buf_offset++;
bufDSend[buf_offset] = r22; buf_offset++;
bufDSend[buf_offset] = r23_a; buf_offset++;

bufDSend[buf_offset] = 0.0; buf_offset++;
bufDSend[buf_offset] = r23_b; buf_offset++;
bufDSend[buf_offset] = r33; buf_offset++;
}

for (iVar = 0; iVar < ICOUNT; iVar++) {
for (iDim = 0; iDim < nDim; iDim++) {
bufDSend[buf_offset] = Cvector(iVar,iDim);
buf_offset++;
}
}

break;

case PERIODIC_LIM_PRIM_1:
case PERIODIC_LIM_SOL_1:

/*--- The first phase of the periodic limiter calculation
ensures that the proper min and max of the solution are found
among all nodes adjacent to periodic faces. ---*/

/*--- We send the min and max over "our" neighbours. ---*/

for (iVar = 0; iVar < ICOUNT; iVar++) {
Sol_Min[iVar] = base_nodes->GetSolution_Min()(iPoint, iVar);
Sol_Max[iVar] = base_nodes->GetSolution_Max()(iPoint, iVar);
}

for (auto jPoint : geometry->nodes->GetPoints(iPoint)) {
for (iVar = 0; iVar < ICOUNT; iVar++) {
Sol_Min[iVar] = min(Sol_Min[iVar], field(jPoint, iVar));
Sol_Max[iVar] = max(Sol_Max[iVar], field(jPoint, iVar));
}
}

for (iVar = 0; iVar < ICOUNT; iVar++) {
bufDSend[buf_offset+iVar] = Sol_Min[iVar];
bufDSend[buf_offset+ICOUNT+iVar] = Sol_Max[iVar];
}

/*--- Rotate the momentum components of the min/max. ---*/

if (rotate_periodic) {
Rotate(zeros, &Sol_Min[1], &bufDSend[buf_offset+1]);
Rotate(zeros, &Sol_Max[1], &bufDSend[buf_offset+ICOUNT+1]);
}

break;

case PERIODIC_LIM_PRIM_2:
case PERIODIC_LIM_SOL_2:

/*--- The second phase of the periodic limiter calculation
ensures that the correct minimum value of the limiter is
found for a node on a periodic face and stores it. ---*/

for (iVar = 0; iVar < ICOUNT; iVar++) {
bufDSend[buf_offset+iVar] = limiter(iPoint, iVar);
}

if (rotate_periodic) {
Rotate(zeros, &limiter(iPoint,1), &bufDSend[buf_offset+1]);
}

break;

default:
SU2_MPI::Error("Unrecognized quantity for periodic communication.",
CURRENT_FUNCTION);
break;
}

Check notice

Code scanning / CodeQL

Long switch case Note

Switch has at least one case that is too long:
PERIODIC_RESIDUAL (68 lines)
Loading
.
Switch has at least one case that is too long:
PERIODIC_LAPLACIAN (58 lines)
Loading
.
Switch has at least one case that is too long:
PERIODIC_SENSOR (50 lines)
Loading
.
Switch has at least one case that is too long:
PERIODIC_PRIM_GG_R (43 lines)
Loading
.
Switch has at least one case that is too long:
PERIODIC_PRIM_ULS_R (150 lines)
Loading
.
Switch has at least one case that is too long:
PERIODIC_LIM_SOL_1 (33 lines)
Loading
.
}
END_SU2_OMP_FOR

Expand Down Expand Up @@ -1233,6 +1255,7 @@
case PERIODIC_SOL_GG_R:
case PERIODIC_PRIM_GG:
case PERIODIC_PRIM_GG_R:
case PERIODIC_AUXVAR_GG:

/*--- For G-G, we accumulate partial gradients then compute
the final value using the entire volume of the periodic cell. ---*/
Expand All @@ -1243,6 +1266,7 @@

break;

case PERIODIC_AUXVAR_LS:
case PERIODIC_SOL_LS: case PERIODIC_SOL_ULS:
case PERIODIC_SOL_LS_R: case PERIODIC_SOL_ULS_R:
case PERIODIC_PRIM_LS: case PERIODIC_PRIM_ULS:
Expand Down Expand Up @@ -2204,19 +2228,19 @@
const auto& solution = base_nodes->GetAuxVar();
auto& gradient = base_nodes->GetAuxVarGradient();

computeGradientsGreenGauss(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_NONE, *geometry,
computeGradientsGreenGauss(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_AUXVAR_GG, *geometry,
*config, solution, 0, base_nodes->GetnAuxVar(), -1, gradient);
}

void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) {
SU2_ZONE_SCOPED

bool weighted = true;
const bool weighted = true;
const auto& solution = base_nodes->GetAuxVar();
auto& gradient = base_nodes->GetAuxVarGradient();
auto& rmatrix = base_nodes->GetRmatrix();

computeGradientsLeastSquares(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config,
computeGradientsLeastSquares(this, MPI_QUANTITIES::AUXVAR_GRADIENT, PERIODIC_AUXVAR_LS, *geometry, *config,
weighted, solution, 0, base_nodes->GetnAuxVar(), -1, gradient, rmatrix);
}

Expand Down
1 change: 1 addition & 0 deletions SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf,
source_cons_jac.resize(nPoint, flamelet_config_options.n_user_scalars) = su2double(0.0);

if (flamelet_config_options.preferential_diffusion) {
nAuxVar = FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS;
AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS) = su2double(0.0);
Grad_AuxVar.resize(nPoint, FLAMELET_PREF_DIFF_SCALARS::N_BETA_TERMS, nDim, 0.0);
}
Expand Down
8 changes: 8 additions & 0 deletions TestCases/hybrid_regression.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -153,6 +153,14 @@ def main():
poiseuille_profile.test_vals_aarch64 = [-12.004276, -7.636719, -0.000000, 2.089953]
test_list.append(poiseuille_profile)

# Streamwise periodic pressure-drop update with two OpenMP threads.
sp_massflow_box = TestCase('sp_massflow_box')
sp_massflow_box.cfg_dir = "incomp_navierstokes/streamwise_periodic"
sp_massflow_box.cfg_file = "massflow_box.cfg"
sp_massflow_box.test_iter = 4
sp_massflow_box.test_vals = [-4.914357, 0.248162, -0.092770]
test_list.append(sp_massflow_box)

# 2D Rotational Periodic
periodic2d = TestCase('periodic2d')
periodic2d.cfg_dir = "navierstokes/periodic2D"
Expand Down
36 changes: 36 additions & 0 deletions TestCases/incomp_navierstokes/streamwise_periodic/massflow_box.cfg
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
% Small generated periodic channel: pressure-drop update must be independent of OpenMP threads.
SOLVER= INC_NAVIER_STOKES
MESH_FORMAT= BOX
MESH_BOX_SIZE= (9, 9, 9)
MESH_BOX_LENGTH= (1, 1, 1)
MESH_BOX_OFFSET= (0, 0, 0)
INC_DENSITY_MODEL= CONSTANT
INC_DENSITY_INIT= 1.0
INC_VELOCITY_INIT= (0.1, 0.0, 0.0)
INC_TEMPERATURE_INIT= 300.0
INC_ENERGY_EQUATION= NO
VISCOSITY_MODEL= CONSTANT_VISCOSITY
MU_CONSTANT= 0.01
MARKER_HEATFLUX= (y_minus,0, y_plus,0, z_minus,0, z_plus,0)
MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)
KIND_STREAMWISE_PERIODIC= MASSFLOW
STREAMWISE_PERIODIC_MASSFLOW= 1.0
STREAMWISE_PERIODIC_PRESSURE_DROP= 0.0
INC_OUTLET_DAMPING= 0.1
NUM_METHOD_GRAD= GREEN_GAUSS
CONV_NUM_METHOD_FLOW= FDS
MUSCL_FLOW= NO
TIME_DISCRE_FLOW= EULER_IMPLICIT
CFL_NUMBER= 1.0
CFL_ADAPT= NO
LINEAR_SOLVER= FGMRES
LINEAR_SOLVER_PREC= ILU
LINEAR_SOLVER_ITER= 10
ITER= 5
CONV_RESIDUAL_MINVAL= -30
OUTPUT_FILES= NONE
OUTPUT_WRT_FREQ= 5
HISTORY_OUTPUT= (RMS_RES, STREAMWISE_PERIODIC)
SCREEN_OUTPUT= (INNER_ITER, RMS_PRESSURE, STREAMWISE_DP, STREAMWISE_MASSFLOW)
SCREEN_WRT_FREQ_INNER= 1
INC_NONDIM= DIMENSIONAL
Loading
Loading