Skip to content

Commit ca3d54e

Browse files
rois1995claudejoshkellyjakbigfooted
authored
Fixes for rotational periodic boundaries: implicit solver, multigrid and limiters (#2961)
## Proposed Changes This is one of **seven related periodic-boundary bug-fix PRs**. | PR | Scope | |---|---| | #2961 (this PR) | Rotational Jacobian, multigrid velocity rotation and limiters | | #2963 | Axis points, slip-wall normals and moving donor volumes | | #2964 | Periodic volumes, neighbor counts, sensors, stencils and coarse flags | | #2965 | Streamwise periodicity, auxiliary gradients and flamelet metadata | | #2966 | Inverse affine transforms in the parser and FEM | | #2968 | Diagnostics for unsupported periodic solver combinations | | #2969 | Translated wall-distance images | The other bug-fix PRs are drafts with their remaining validation and implementation limits documented. Their grouping does not imply a linear merge order; #2963 depends on this PR. The fully coupled implicit operator proposed in #2967 is a separate experimental enhancement. It has been closed and deferred because no overall benefit has been established and convergence/validation issues remain. This PR retains the existing implicit approximation. Three fixes for rotational periodic boundaries, one commit each. The test cases are in the comment below. - **Implicit runs:** the Jacobian block that a periodic point gets from its match was rotated on one side only. It is now rotated on both sides. - **Multigrid:** the coarse levels did not rotate the velocity, so the run did not converge. The velocity is now rotated on all levels. - **Limiters:** the limiters of the velocity components were rotated like a vector, but they are separate numbers. On the periodic boundary one velocity component became first order, and limiters could be negative. Each side now sends the velocity bounds and reconstruction increments in the frame of the other side, so each limiter uses the complete stencil. Nothing changes for a translation. New tests `periodic2d_no_limiter` and `periodic2d_multigrid`; `periodic2d` is now also in the MPI regression. Their x86 reference values pass in CI. Existing reference values that change: `aachen_turbine_restart`, `jones_turbocharger_restart` and `multi_interface` (rotational periodic, the Jacobian block), and the hybrid `periodic2d` (the Jacobian block and the limiters). The x86 periodic tests and the other changed cases pass in CI; the Jones restart residual references were refreshed from the serial, MPI and hybrid CI logs. The aarch64 references still need the CI runs. No conflict with #2945 and #2959. #2956 also changes the reference values of `jones_turbocharger_restart`; whichever is merged second takes them from the CI. ## Related Work - #652 added the present implementation of the periodic boundaries. - #1563 fixed the rotation in 2D and added the `periodic2D` case used here. - #1467, #763 and #1585 discuss the incomplete linear system at periodic points; this PR keeps the existing approximation. ## PR Checklist - [x] I am submitting my contribution to the develop branch. - [x] My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson). - [x] My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/). - [x] I used the pre-commit hook to prevent dirty commits and used `pre-commit run --all` to format old commits. - [x] I have added a test case that demonstrates my contribution, if necessary. - [ ] I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary. --------- Co-authored-by: rois1995 <rauxa1995@hotmail.it> Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> Co-authored-by: Josh Kelly <81244680+joshkellyjak@users.noreply.github.com> Co-authored-by: Nijso <bigfootedrockmidget@hotmail.com>
1 parent 7c18ea5 commit ca3d54e

15 files changed

Lines changed: 621 additions & 31 deletions

File tree

‎Common/include/toolboxes/geometry_toolbox.hpp‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
#pragma once
2828

2929
#include <cmath>
30+
#include <algorithm>
3031

3132
namespace GeometryToolbox {
3233
/// \addtogroup GeometryToolbox
@@ -215,6 +216,31 @@ inline void Rotate(const Scalar R[][nDim], const Scalar* O, const Scalar* d, Sca
215216
}
216217
}
217218

219+
/*! \return Whether any of the three supplied rotation angles is nonzero. */
220+
template <class Scalar>
221+
inline bool HasRotation(const Scalar* angles) {
222+
// Configured zero stays exact after degree conversion and negation; retain every nonzero rotation.
223+
return angles[0] != 0.0 || angles[1] != 0.0 || angles[2] != 0.0;
224+
}
225+
226+
/*! \brief Rotate component bounds in place, enclosing the rotated box. */
227+
template <class Scalar, int nDim>
228+
inline void RotateBox(const Scalar R[][nDim], Scalar* vMin, Scalar* vMax) {
229+
using std::max;
230+
using std::min;
231+
Scalar rotMin[nDim] = {0.0}, rotMax[nDim] = {0.0};
232+
for (int iDim = 0; iDim < nDim; ++iDim) {
233+
for (int jDim = 0; jDim < nDim; ++jDim) {
234+
const Scalar fromMin = R[iDim][jDim] * vMin[jDim];
235+
const Scalar fromMax = R[iDim][jDim] * vMax[jDim];
236+
rotMin[iDim] += min(fromMin, fromMax);
237+
rotMax[iDim] += max(fromMin, fromMax);
238+
}
239+
}
240+
std::copy_n(rotMin, nDim, vMin);
241+
std::copy_n(rotMax, nDim, vMax);
242+
}
243+
218244
/*! \brief Tangent projection */
219245
template <class Mat, class Scalar, class Int>
220246
inline void TangentProjection(Int nDim, const Mat& tensor, const Scalar* vector, Scalar* proj) {

‎SU2_CFD/include/limiters/CLimiterDetails.hpp‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,16 @@ struct LimiterHelpers
7272
{
7373
FORCEINLINE static Type epsilon() {return std::numeric_limits<passivedouble>::epsilon();}
7474

75+
/*! \brief MUSCL reconstruction increment using the displacement to the middle of an edge. */
76+
template <class Int>
77+
FORCEINLINE static Type reconstructionIncrement(Int nDim, const Type* halfEdge, const Type* gradient,
78+
const Type& value_i, const Type& value_j, const Type& kappa) {
79+
Type proj = 0.0;
80+
for (Int iDim = 0; iDim < nDim; ++iDim) proj += halfEdge[iDim] * gradient[iDim];
81+
const Type cent = 0.5 * (value_j - value_i);
82+
return umusclProjection(proj, cent, kappa);
83+
}
84+
7585
FORCEINLINE static Type umusclProjection(const Type& grad_proj, const Type& delta, const Type& kappa)
7686
{
7787
/*-------------------------------------------------------------------*/

‎SU2_CFD/include/limiters/computeLimiters_impl.hpp‎

Lines changed: 40 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -109,11 +109,29 @@ void computeLimiters_impl(CSolver* solver,
109109

110110
limiterDetails.preprocess(geometry, config, varBegin, varEnd, field);
111111

112+
/*--- With rotational periodicity, the first periodic comm. also brings the min/max
113+
* projections over the edges of the periodic matches (stored after each other),
114+
* because the limiters of the velocity cannot be compared across a rotation. ---*/
115+
116+
su2activematrix* periodicProj = nullptr;
117+
118+
if (periodic && (kindPeriodicComm1 == PERIODIC_LIM_PRIM_1))
119+
periodicProj = solver->GetPeriodicProjections(geometry, config);
120+
112121
/*--- Initialize all min/max field values if we have
113122
* periodic comms. otherwise do it inside main loop. ---*/
114123

115124
if (periodic)
116125
{
126+
if (periodicProj != nullptr)
127+
{
128+
SU2_OMP_FOR_STAT(chunkSize)
129+
for (auto iPoint = 0ul; iPoint < periodicProj->rows(); ++iPoint)
130+
for (auto iVar = 0ul; iVar < periodicProj->cols(); ++iVar)
131+
(*periodicProj)(iPoint,iVar) = 0.0;
132+
END_SU2_OMP_FOR
133+
}
134+
117135
SU2_OMP_FOR_STAT(chunkSize)
118136
for (size_t iPoint = 0; iPoint < nPoint; ++iPoint)
119137
for (size_t iVar = varBegin; iVar < varEnd; ++iVar)
@@ -165,6 +183,23 @@ void computeLimiters_impl(CSolver* solver,
165183
for (size_t iVar = varBegin; iVar < varEnd; ++iVar)
166184
projMax[iVar] = projMin[iVar] = 0.0;
167185

186+
if (periodicProj != nullptr && nodes->GetPeriodicBoundary(iPoint))
187+
{
188+
/*--- Start from the min/max over the edges of the periodic matches. ---*/
189+
190+
const auto* projections = solver->GetPeriodicProjection(iPoint);
191+
if (projections != nullptr) {
192+
for (auto iVar = varBegin; iVar < varEnd; ++iVar) {
193+
const auto& periodicMin = projections[iVar];
194+
const auto& periodicMax = projections[periodicProj->cols()/2 + iVar];
195+
AD::SetPreaccIn(periodicMin);
196+
AD::SetPreaccIn(periodicMax);
197+
projMin[iVar] = periodicMin;
198+
projMax[iVar] = periodicMax;
199+
}
200+
}
201+
}
202+
168203
/*--- Compute max/min projection and values over direct neighbors. ---*/
169204

170205
for (auto jPoint : geometry.nodes->GetPoints(iPoint)) {
@@ -175,22 +210,16 @@ void computeLimiters_impl(CSolver* solver,
175210
/*--- Distance vector from iPoint to face (middle of the edge). ---*/
176211

177212
su2double dist_ij[nDim] = {0.0};
178-
179-
for(size_t iDim = 0; iDim < nDim; ++iDim)
213+
for (size_t iDim = 0; iDim < nDim; ++iDim)
180214
dist_ij[iDim] = 0.5 * (coord_j[iDim] - coord_i[iDim]);
181215

182216
/*--- Project each variable, update min/max. ---*/
183217

184218
for(size_t iVar = varBegin; iVar < varEnd; ++iVar)
185219
{
186-
su2double proj = 0.0;
187-
188-
for(size_t iDim = 0; iDim < nDim; ++iDim)
189-
proj += dist_ij[iDim] * gradient(iPoint,iVar,iDim);
190-
191220
AD::SetPreaccIn(field(jPoint,iVar));
192-
const su2double cent = 0.5 * (field(jPoint,iVar) - field(iPoint,iVar));
193-
proj = LimiterHelpers<>::umusclProjection(proj, cent, umusclKappa);
221+
const su2double proj = LimiterHelpers<>::reconstructionIncrement(nDim, dist_ij,
222+
gradient[iPoint][iVar], field(iPoint,iVar), field(jPoint,iVar), umusclKappa);
194223

195224
projMax[iVar] = max(projMax[iVar], proj);
196225
projMin[iVar] = min(projMin[iVar], proj);
@@ -224,7 +253,8 @@ void computeLimiters_impl(CSolver* solver,
224253
}
225254
END_SU2_OMP_FOR
226255

227-
/*--- Account for periodic effects, take the minimum limiter on each periodic pair. ---*/
256+
/*--- Account for periodic effects, take the minimum limiter on each periodic pair
257+
* (except for the velocity with rotational periodicity). ---*/
228258
if (periodic)
229259
{
230260
for (size_t iPeriodic = 1; iPeriodic <= config.GetnMarker_Periodic()/2; ++iPeriodic)

‎SU2_CFD/include/solvers/CFVMFlowSolverBase.inl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -261,7 +261,7 @@ void CFVMFlowSolverBase<V, R>::CommunicateInitialState(CGeometry* geometry, cons
261261
CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS);
262262
}
263263
SetImplicitPeriodic(euler_implicit);
264-
if (MGLevel == MESH_0) SetRotatePeriodic(true);
264+
SetRotatePeriodic(true);
265265

266266
/*--- Perform the MPI communication of the solution ---*/
267267

‎SU2_CFD/include/solvers/CSolver.hpp‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@
3737
#include <algorithm>
3838
#include <iostream>
3939
#include <set>
40+
#include <unordered_map>
4041
#include <stdlib.h>
4142
#include <stdio.h>
4243

@@ -146,6 +147,8 @@ class CSolver {
146147

147148
bool rotate_periodic; /*!< \brief Flag that controls whether the periodic solution needs to be rotated for the solver. */
148149
bool implicit_periodic; /*!< \brief Flag that controls whether the implicit system should be treated by the periodic BC comms. */
150+
su2activematrix PeriodicProj; /*!< \brief Min/max reconstruction increments at periodic receive points (for limiters). */
151+
std::unordered_map<unsigned long, size_t> PeriodicProjIndex; /*!< \brief Local point to compact projection row. */
149152

150153
bool dynamic_grid; /*!< \brief Flag that determines whether the grid is dynamic (moving or deforming + grid velocities). */
151154

@@ -232,6 +235,8 @@ class CSolver {
232235
return base_nodes;
233236
}
234237

238+
bool HasPeriodicRotation(const CConfig& config) const;
239+
235240
/*!
236241
* \brief Helper function to define the type and number of variables per point for each communication type.
237242
* \param[in] config - Definition of the particular problem.
@@ -4231,6 +4236,22 @@ class CSolver {
42314236
*/
42324237
inline void SetRotatePeriodic(bool val_rotate_periodic) { rotate_periodic = val_rotate_periodic; }
42334238

4239+
/*!
4240+
* \brief Storage for the limiters with rotational periodicity: the min and max, over the edges of the periodic
4241+
* matches of each point, of the reconstruction increments (communicated with PERIODIC_LIM_PRIM_1).
4242+
* \param[in] geometry - Periodic receive points of this mesh level.
4243+
* \param[in] config - Definition of the particular problem.
4244+
* \return The matrix (unique periodic receive points x 2*nPrimVarGrad, min then max), nullptr without rotation.
4245+
*/
4246+
su2activematrix* GetPeriodicProjections(const CGeometry& geometry, const CConfig& config);
4247+
4248+
/*! \brief Reconstruction increment bounds for a periodic receive point, nullptr for other points. */
4249+
inline su2double* GetPeriodicProjection(unsigned long iPoint) {
4250+
const auto& indices = PeriodicProjIndex;
4251+
const auto it = indices.find(iPoint);
4252+
return it == indices.end() ? nullptr : PeriodicProj[it->second];
4253+
}
4254+
42344255
/*!
42354256
* \brief Retrieve the solver name for output purposes.
42364257
* \returns Name of the solver.

0 commit comments

Comments
 (0)