Skip to content
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.

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

Correction of symmetry plane implementation #2194

Open
wants to merge 177 commits into
base: develop
Choose a base branch
from

Conversation

bigfooted
Copy link
Contributor

@bigfooted bigfooted commented Jan 3, 2024

Proposed Changes

The implementation of the symmetry plane is incomplete. We follow here the book of Blazek, Computational Fluid Dynamics: Principles and Applications . According to Blazek, (chapter 8.6) 4 conditions have to be met on a symmetry plane:

  1. No flux across the boundary
  2. $\overline n \cdot \nabla \phi = 0$ (gradient of scalars are zero)
  3. $\overline n \cdot \nabla (\overline U \cdot \overline t) = 0$ (gradient of tangential velocity normal to boundary)
  4. $\overline t \cdot \nabla (\overline U \cdot \overline n) = 0$ (gradient of normal velocity along the boundary)

Points 2-4 all deal with gradients and can be dealt with in the gradient computation, i.e. Green-Gauss or Least Squares. According to Blazek, 2 approaches can be followed. "One possibility is to construct the missing half of the control volume by mirroring the grid on the boundary. The fluxes and gradients are then evaluated like in the interior using reflected flow variables." This approach can be implemented in an easy way when computing the Green-Gauss gradients. In SU2, routines are already in place that deal with GG gradients on boundaries. Here, we just have to identify the symmetry planes and mirror the flux through the faces.
Blazek continues: "The second methodology computes the fluxes for the halved control volume (but not accross the boundary). The components of the residual normal to the symmetry plane are then zeroed out. It is also necessary to correct normal vectors of those faces of the control volume, which touch the boundary. The modification consists of removing all components of the face vector, which are normal to the symmetry plane. The gradients also have to be corrected according to Eq. (8.40) [points 2,3,4 above]"

  1. Modify Green-Gauss and Least-Squares gradient computation for symmetry planes
    -> This seems to work fine
  2. Fix normal components of edges connected to symmetry plane
    -> This seems to work fine
  3. Fix residual components of velocity on the symmetry plane
    -> This seems to work fine
  4. Multiple symmetry planes, including Euler walls
    -> This seems to work fine
  5. multigrid
    -> Works fine now that we use the agglomeration as in the paper of Diskin.
  6. grid movement
  7. -> Works fine with the new fix to the residuals (energy update)

Related Work

#1168

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

  • I am submitting my contribution to the develop branch.
  • My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
  • My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
  • I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
  • 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.

@bigfooted bigfooted marked this pull request as draft January 3, 2024 23:34
Common/src/geometry/CPhysicalGeometry.cpp Fixed Show fixed Hide fixed
Common/src/geometry/CPhysicalGeometry.cpp Fixed Show fixed Hide fixed
Common/src/geometry/CPhysicalGeometry.cpp Fixed Show fixed Hide fixed
Common/src/geometry/CPhysicalGeometry.cpp Fixed Show fixed Hide fixed
Common/src/geometry/CPhysicalGeometry.cpp Fixed Show fixed Hide fixed
Common/include/CConfig.hpp Outdated Show resolved Hide resolved
SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp Outdated Show resolved Hide resolved
SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp Outdated Show resolved Hide resolved
SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp Outdated Show resolved Hide resolved
SU2_CFD/include/gradients/computeGradientsSymmetry.hpp Outdated Show resolved Hide resolved
SU2_CFD/include/gradients/computeGradientsSymmetry.hpp Outdated Show resolved Hide resolved
Comment on lines +70 to +73
for (auto iDim = 0u; iDim < nDim; iDim++) {
for (auto jDim = 0u; jDim < nDim; jDim++) {
for (auto kDim = 0u; kDim < nDim; kDim++) {
for (auto mDim = 0u; mDim < nDim; mDim++) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

4 loops may not be the most efficient way of doing 2 matrix multiplications.
From the way the basis/tensor is constructed there are probably some nice cancelations that make the correction simpler.

Comment on lines +295 to +301
su2activematrix Gradients_iPoint(varEnd-varBegin,nDim);

for (auto iVar = varBegin; iVar < varEnd; iVar++) {
for (auto iDim = 0u; iDim < nDim; iDim++) {
Gradients_iPoint[iVar][iDim] = gradient(iPoint, iVar, iDim);
}
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can probably use the gradient iterator to avoid the copy.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What do you mean with the gradient iterator?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The 3D matrix returns an object that looks like a matrix but does not require copying the data. I can help with that sort of thing once all the tests are green.

Comment on lines +128 to +129
/*--- map transpose T' * grad(phi) ---*/
gradPhiReflected[jDim] += TensorMap[jDim][iDim] * Gradients_iPoint[iVar][iDim];
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not the map transposed right? You are adding over iDim

Comment on lines +135 to +143
/*--- gradient in direction normal to symmetry is cancelled ---*/
gradPhiReflected[0] = 0.0;

/*--- Now transform back ---*/
for (auto jDim = 0u; jDim < nDim; jDim++) {
for (auto iDim = 0u; iDim < nDim; iDim++) {
gradPhi[jDim] += TensorMap[iDim][jDim] * gradPhiReflected[iDim];
}
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For example, for scalars these operations are equivalent to:

g = T' * [0 0 0; 0 1 0; 0 0 1] * T * g

(setting the first component to 0 in matrix form).
This is equivalent to

g = (I - T(1,:)' * T(1,:)) * g

(which is the same as taking out the normal component)
The tangential and orthogonal directions cancel out (as desired since they are arbitrary).

After you review the tests we can look into these simplifications.

Comment on lines +1151 to +1159
// unsigned short Syms[MAXNSYMS] = {0};
// unsigned short nSym = 0;
// for (size_t iMarker = 0; iMarker < geometry->GetnMarker(); ++iMarker) {
// if ((config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE) ||
// (config->GetMarker_All_KindBC(iMarker) == EULER_WALL)) {
// Syms[nSym] = iMarker;
// nSym++;
// }
// }

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
Comment on lines +1207 to +1234
// if (iPoint>1500) {
// cout << "normal of the primary: "<<UnitNormalPrim[0] << " " << UnitNormalPrim[1] << " " << UnitNormalPrim[2] << endl;
// cout << "normal of the secondary:"<<UnitNormal[0] << " " << UnitNormal[1] << " " << UnitNormal[2] << endl;
// }
// /*--- Correct the current normal as n2_new = n2 - (n2.n1)n1 ---*/
// su2double ProjNorm = 0.0;
// for (auto iDim = 0u; iDim < nDim; iDim++) ProjNorm += UnitNormal[iDim] * UnitNormalPrim[iDim];
// /*--- We check if the normal of the 2 planes coincide.
// * We only update the normal if the normals of the symmetry planes are different. ---*/
// if (fabs(1.0-ProjNorm)>EPS) {
// for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] -= ProjNorm * UnitNormalPrim[iDim];
// /*--- Make normalized vector ---*/
// su2double newarea=GeometryToolbox::Norm(nDim, UnitNormal);
// for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = UnitNormal[iDim]/newarea;
// if (iPoint>1500)
// cout << " ipoint=" <<iPoint << " " << val_marker << " " << Syms[iMarker] << " " << geometry->nodes->GetCoord(iPoint, 0) << " " <<
// geometry->nodes->GetCoord(iPoint, 1) << " " <<
// geometry->nodes->GetCoord(iPoint, 2) << " , n="
// << UnitNormal[0] << " "
// << UnitNormal[1] << " "
// << UnitNormal[2] << endl;
// }
// }
// }
// }
// }
// }
// }

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
Comment on lines +1202 to +1205
// su2double AreaPrim = GeometryToolbox::Norm(nDim, NormalPrim);
// for(unsigned short iDim = 0; iDim < nDim; iDim++) {
// UnitNormalPrim[iDim] = NormalPrim[iDim] / AreaPrim;
// }

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
Comment on lines +1190 to +1195
// /*--- Loop over all points on the other symmetry and check if current iPoint is on the symmetry ---*/
// for (auto jVertex = 0ul; jVertex < geometry->nVertex[Syms[iMarker]]; jVertex++) {
// const auto jPoint = geometry->vertex[Syms[iMarker]][jVertex]->GetNode();
// if (!geometry->nodes->GetDomain(jPoint)) continue;
// if (iPoint==jPoint) {
// //cout << "marker = " << val_marker<< ", other Marker =" << Syms[iMarker] << endl;

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
Comment on lines +1181 to +1187
// if (nSym>1) {
// /*--- Normal of the primary symmetry plane ---*/
// su2double NormalPrim[MAXNDIM] = {0.0}, UnitNormalPrim[MAXNDIM] = {0.0};
// /*--- Step 2: are we on a shared node? ---*/
// for (auto iMarker=0;iMarker<nSym;iMarker++) {
// /*--- We do not want the current symmetry ---*/
// if (val_marker!= Syms[iMarker]) {

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
/*--- Allocation of variables necessary for viscous fluxes. ---*/
su2double ProjGradient, ProjNormVelGrad, ProjTangVelGrad, TangentialNorm,
Tangential[MAXNDIM] = {0.0}, GradNormVel[MAXNDIM] = {0.0}, GradTangVel[MAXNDIM] = {0.0};
// static constexpr size_t MAXNSYMS = 100;

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
3 participants