Skip to content

Commit

Permalink
fixes to face normal correction on symmetry
Browse files Browse the repository at this point in the history
  • Loading branch information
@bigfooted
bigfooted committed Jan 10, 2024
1 parent a3e7221 commit 48f6d3d
Show file tree
Hide file tree
Showing 5 changed files with 282 additions and 277 deletions.
106 changes: 52 additions & 54 deletions Common/src/geometry/CPhysicalGeometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7282,60 +7282,58 @@ void CPhysicalGeometry::SetBoundControlVolume(const CConfig* config, unsigned sh
* 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. ---*/
// SU2_OMP_FOR_DYN(1)
// for (unsigned short iMarker = 0; iMarker < nMarker; iMarker++) {
// if (config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE) {
// for (unsigned long iElem = 0; iElem < nElem_Bound[iMarker]; iElem++) {

// const auto nNodes = bound[iMarker][iElem]->GetnNodes();
// for (unsigned short iNode = 0; iNode < nNodes; iNode++) {

// const auto iPoint = bound[iMarker][iElem]->GetNode(iNode);
// const auto iVertex = nodes->GetVertex(iPoint, iMarker);

// const auto Normal_Sym = vertex[iMarker][iVertex]->GetNormal();
// su2double Area = GeometryToolbox::Norm(nDim, Normal_Sym);
// su2double UnitNormal_Sym[MAXNDIM] = {0.0};

// for(unsigned short iDim = 0; iDim<nDim; iDim++){
// UnitNormal_Sym[iDim] = Normal_Sym[iDim]/Area;
// }

// for (unsigned short iNeigh = 0; iNeigh < bound[iMarker][iElem]->GetnNeighbor_Nodes(iNode); ++iNeigh){
// su2double Product = 0.0;
// unsigned long jPoint = nodes->GetPoint(iPoint,iNeigh);
// /*---Check if neighbour point is on the same plane as the symmetry plane
// by computing the internal product and of the Normal Vertex vector and
// the vector connecting iPoint and jPoint. If the product is lower than
// estabilished tolerance (to account for Numerical errors) both points are
// in the same plane as SYMMETRY_PLANE---*/
// su2double Tangent[MAXNDIM]={0.0};
// for(unsigned short iDim = 0; iDim<nDim; iDim++){
// Tangent[iDim] = nodes->GetCoord(jPoint,iDim) - nodes->GetCoord(iPoint,iDim);
// Product += Tangent[iDim] * Normal_Sym[iDim];
// }

// if (abs(Product) < EPS) {
// Product = 0.0;

// unsigned long iEdge = nodes->GetEdge(iPoint,iNeigh);
// su2double Normal[MAXNDIM] = {0.0};
// edges->GetNormal(iEdge,Normal);

// for(unsigned short iDim = 0; iDim<nDim; iDim++)
// Product += Normal[iDim]*UnitNormal_Sym[iDim];

// for(unsigned short iDim = 0; iDim<nDim; iDim++)
// Normal[iDim]-=Product*UnitNormal_Sym[iDim];

// edges->SetNormal(iEdge,Normal);
// }
// }
// }
// } // loop over elements
// } // if symmetry
// } // loop over markers
// END_SU2_OMP_FOR
SU2_OMP_FOR_DYN(1)
for (unsigned short iMarker = 0; iMarker < nMarker; iMarker++) {
cout << "loop over markers" << endl;
if (config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE) {
for (unsigned long iVertex = 0; iVertex < GetnVertex(iMarker); iVertex++) {

const auto iPoint = vertex[iMarker][iVertex]->GetNode();

const auto Normal_Sym = vertex[iMarker][iVertex]->GetNormal();

su2double Area = GeometryToolbox::Norm(nDim, Normal_Sym);
su2double UnitNormal_Sym[MAXNDIM] = {0.0};

for (unsigned short iDim = 0; iDim < nDim; iDim++) {
UnitNormal_Sym[iDim] = Normal_Sym[iDim]/Area;
}

for (unsigned short iNeigh = 0; iNeigh < nodes->GetnPoint(iPoint); ++iNeigh){
su2double Product = 0.0;
unsigned long jPoint = nodes->GetPoint(iPoint,iNeigh);
/*---Check if neighbour point is on the same plane as the symmetry plane
by computing the internal product of the Normal Vertex vector and
the vector connecting iPoint and jPoint. If the product is lower than
estabilished tolerance (to account for Numerical errors) both points are
in the same plane as SYMMETRY_PLANE---*/
su2double Tangent[MAXNDIM]={0.0};
for(unsigned short iDim = 0; iDim<nDim; iDim++){
Tangent[iDim] = nodes->GetCoord(jPoint,iDim) - nodes->GetCoord(iPoint,iDim);
Product += Tangent[iDim] * Normal_Sym[iDim];
}

if (abs(Product) < 1.0e-8) {
cout << "iPoint " << iPoint << " and jPoint " << jPoint << " are in the symmetry plane " << endl;
Product = 0.0;

unsigned long iEdge = nodes->GetEdge(iPoint,iNeigh);
su2double Normal[MAXNDIM] = {0.0};
edges->GetNormal(iEdge,Normal);
cout << "normal="<<Normal[0] << " " << Normal[1] << endl;
for(unsigned short iDim = 0; iDim<nDim; iDim++)
Product += Normal[iDim]*UnitNormal_Sym[iDim];

for(unsigned short iDim = 0; iDim<nDim; iDim++)
Normal[iDim]-=Product*UnitNormal_Sym[iDim];

edges->SetNormal(iEdge,Normal);
} // if in-plane of symmetry
} // loop over neighbors
} // loop over vertices
} // if symmetry
} // loop over markers
END_SU2_OMP_FOR
}

void CPhysicalGeometry::VisualizeControlVolume(const CConfig* config) const {
Expand Down
14 changes: 1 addition & 13 deletions SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -224,6 +224,7 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
(config.GetMarker_All_KindBC(iMarker) != NEARFIELD_BOUNDARY) &&
(config.GetMarker_All_KindBC(iMarker) != SYMMETRY_PLANE) &&
(config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) {

/*--- Work is shared in inner loop as two markers
* may try to update the same point. ---*/

Expand All @@ -239,18 +240,9 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
su2double volume = nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint);
const auto area = geometry.vertex[iMarker][iVertex]->GetNormal();


// When the node is shared with a symmetry we need to mirror the contribution of
// the face that is coincident with the inlet/outlet
if (nodes->GetSymmetry(iPoint) && nodes->Getinoutfar(iPoint)) {
cout << "iPoint "
<< iPoint
<< " is on a symmetry plane and an inlet/outlet"
<< nodes->GetSymmetry(iPoint)
<< ", "
<< nodes->Getinoutfar(iPoint)<< endl;

cout << " face area = " << area[0] <<", " << area[1] << endl;

// we have to find the edges that were missing in the symmetry computations.
// So we find the jPoints that are on the inlet plane
Expand Down Expand Up @@ -278,7 +270,6 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
cout << "flux ("<<iPoint<<","<<iVar<<")="<<flux[iVar]<<endl;
}


// find the point iPoint on the symmetry plane and get the vertex normal at ipoint wrt the symmetry plane
const su2double* VertexNormal = getVertexNormalfromPoint(config, geometry,iPoint);
cout << " vertex normal = " << VertexNormal[0] <<", " << VertexNormal[1] << endl;
Expand Down Expand Up @@ -309,9 +300,6 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
fluxReflected[iDim + 1] = flux[iDim + 1] - 2.0 * ProjFlux * UnitNormal[iDim];
}




for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
cout << "fluxReflected ("<<iVar<<")="<<fluxReflected[iVar]<<endl;

Expand Down
Loading

0 comments on commit 48f6d3d

Please sign in to comment.