Skip to content

Commit

Permalink
corrected GG gradient computation on shared nodes
Browse files Browse the repository at this point in the history
  • Loading branch information
@bigfooted
bigfooted committed Jan 8, 2024
1 parent 2840592 commit 7387f8c
Show file tree
Hide file tree
Showing 6 changed files with 206 additions and 65 deletions.
7 changes: 7 additions & 0 deletions Common/include/CConfig.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6403,6 +6403,13 @@ class CConfig {
*/
bool Getinoutfar(unsigned short iMarker) const;

/*!
* \brief Determines whether a marker with index iMarker is a symmetry.
* \param iMarker
* \return <TRUE> it marker with index iMarker is a solid boundary.
*/
bool GetSymmetry(unsigned short iMarker) const;

/*!
* \brief Determines whether a marker with index iMarker is a viscous no-slip boundary.
* \param iMarker
Expand Down
3 changes: 3 additions & 0 deletions Common/include/geometry/dual_grid/CPoint.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ class CPoint {
su2vector<bool>
SolidBoundary; /*!< \brief To see if a point belong to the physical boundary (without includin MPI). */
su2vector<bool> inoutfar;
su2vector<bool> Symmetry;
su2vector<bool>
ViscousBoundary; /*!< \brief To see if a point belong to the physical boundary (without includin MPI). */
su2vector<bool>
Expand Down Expand Up @@ -385,7 +386,9 @@ class CPoint {

// nijso: temporary
inline void Setinoutfar(unsigned long iPoint, bool boundary) { inoutfar(iPoint) = boundary; }
inline void SetSymmetry(unsigned long iPoint, bool boundary) { Symmetry(iPoint) = boundary; }
inline bool Getinoutfar(unsigned long iPoint) const { return inoutfar(iPoint); }
inline bool GetSymmetry(unsigned long iPoint) const { return Symmetry(iPoint); }

/*!
* \brief Set if a point belong to the boundary.
Expand Down
8 changes: 4 additions & 4 deletions Common/src/CConfig.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7936,10 +7936,10 @@ bool CConfig::GetSolid_Wall(unsigned short iMarker) const {

// nijso: temporary
bool CConfig::Getinoutfar(unsigned short iMarker) const {

return (Marker_All_KindBC[iMarker] == INLET_FLOW ||
Marker_All_KindBC[iMarker] == OUTLET_FLOW ||
Marker_All_KindBC[iMarker] == FAR_FIELD);
return (Marker_All_KindBC[iMarker] == INLET_FLOW || Marker_All_KindBC[iMarker] == OUTLET_FLOW);
}
bool CConfig::GetSymmetry(unsigned short iMarker) const {
return (Marker_All_KindBC[iMarker] == SYMMETRY_PLANE);
}

void CConfig::SetSurface_Movement(unsigned short iMarker, unsigned short kind_movement) {
Expand Down
8 changes: 7 additions & 1 deletion Common/src/geometry/CPhysicalGeometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3422,7 +3422,10 @@ void CPhysicalGeometry::SetBoundaries(CConfig* config) {

if (config->GetViscous_Wall(iMarker)) nodes->SetViscousBoundary(Point_Surface, true);
// nijso: temporary
if (config->Getinoutfar(iMarker)) nodes->Setinoutfar(Point_Surface, true);
// cout << "nijso: setting inoutfar for ipoint "<<Point_Surface <<" on marker "<<
// config->GetMarker_All_KindBC(iMarker) << endl; if (config->Getinoutfar(iMarker)) {cout << "inoutfar
// found."<<endl;nodes->Setinoutfar(Point_Surface, true);} if (config->GetSymmetry(iMarker)) {cout << "symmetry
// found"<<endl;nodes->SetSymmetry(Point_Surface, true);}

if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) nodes->SetPeriodicBoundary(Point_Surface, true);
}
Expand Down Expand Up @@ -4601,6 +4604,9 @@ void CPhysicalGeometry::SetRCM_Ordering(CConfig* config) {
if (config->GetSolid_Wall(iMarker)) nodes->SetSolidBoundary(InvResult[iPoint], true);

if (config->GetViscous_Wall(iMarker)) nodes->SetViscousBoundary(InvResult[iPoint], true);
// nijso: temporary
if (config->Getinoutfar(iMarker)) nodes->Setinoutfar(InvResult[iPoint], true);
if (config->GetSymmetry(iMarker)) nodes->SetSymmetry(InvResult[iPoint], true);

if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY)
nodes->SetPeriodicBoundary(InvResult[iPoint], true);
Expand Down
1 change: 1 addition & 0 deletions Common/src/geometry/dual_grid/CPoint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) {
Boundary.resize(npoint) = false;
SolidBoundary.resize(npoint) = false;
inoutfar.resize(npoint) = false; // nijso temporary
Symmetry.resize(npoint) = false; // nijso temporary
ViscousBoundary.resize(npoint) = false;
PhysicalBoundary.resize(npoint) = false;
PeriodicBoundary.resize(npoint) = false;
Expand Down
244 changes: 184 additions & 60 deletions SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,22 @@

namespace detail {

// find local vertex on a symmetry marker using global iPoint
inline su2double* getVertexNormalfromPoint(const CConfig& config, CGeometry& geometry, unsigned long iPointGlobal){
unsigned long iPointSym=0;
for (size_t iMarker = 0; iMarker < geometry.GetnMarker(); ++iMarker) {
if (config.GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE) {
for (size_t iVertex = 0; iVertex < geometry.GetnVertex(iMarker); ++iVertex) {
iPointSym = geometry.vertex[iMarker][iVertex]->GetNode();
if (iPointSym == iPointGlobal)
return geometry.vertex[iMarker][iVertex]->GetNormal();
}
}
}
cout << "point is not found " << endl;
exit(0);
}

/*!
* \brief Compute the gradient of a field using the Green-Gauss theorem.
* \ingroup FvmAlgos
Expand Down Expand Up @@ -144,13 +160,13 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
for (size_t iNeigh = 0; iNeigh < nodes->GetnPoint(iPoint); ++iNeigh) {
size_t iEdge = nodes->GetEdge(iPoint, iNeigh);
size_t jPoint = nodes->GetPoint(iPoint, iNeigh);
su2double* icoord = nodes->GetCoord(iPoint);
su2double* jcoord = nodes->GetCoord(jPoint);
su2double* icoord = nodes->GetCoord(iPoint); // nijso: debugging

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable icoord is not used.
su2double* jcoord = nodes->GetCoord(jPoint); // nijso: debugging

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable jcoord is not used.

su2double Veli[nDim] = {0.0};
for (size_t iDim = 0; iDim < nDim; ++iDim) Veli[iDim] = field(iPoint, iDim + 1);
su2double Velj[nDim] = {0.0};
for (size_t iDim = 0; iDim < nDim; ++iDim) Velj[iDim] = field(jPoint, iDim + 1);
//su2double Veli[nDim] = {0.0};
//for (size_t iDim = 0; iDim < nDim; ++iDim) Veli[iDim] = field(iPoint, iDim + 1);
//su2double Velj[nDim] = {0.0};
//for (size_t iDim = 0; iDim < nDim; ++iDim) Velj[iDim] = field(jPoint, iDim + 1);

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

/*--- Determine if edge points inwards or outwards of iPoint.
* If inwards we need to flip the area vector. ---*/
Expand Down Expand Up @@ -191,53 +207,88 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
}
}

/*--- For nodes shared with walls, we can simply add the mirrored contribution. The nonmirrored
* contributions is added
* in the routine below. ---*/
if (nodes->GetSolidBoundary(iPoint)) {
su2double volume = nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint);
/*--- First, use the values at node i only (better to use entire face but we do not have it) ---*/
for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
flux[iVar] = field(iPoint, iVar) / volume;
fluxReflected[iVar] = flux[iVar];
}
/*--- project the flux ---*/
su2double ProjFlux = 0.0;
for (size_t iDim = 0; iDim < nDim; iDim++) ProjFlux += flux[iDim + 1] * UnitNormal[iDim];

for (size_t iDim = 0; iDim < nDim; iDim++)
fluxReflected[iDim + 1] = flux[iDim + 1] - 2.0 * ProjFlux * UnitNormal[iDim];

for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
for (size_t iDim = 0; iDim < nDim; ++iDim) {
gradient(iPoint, iVar, iDim) -= fluxReflected[iVar] * areaReflected[iDim];
}
}
}

/*--- check if point is shared with an inlet, outlet or far_field ---*/
if (nodes->Getinoutfar(iPoint)) {
su2double volume = nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint);
/*--- First, use the values at node i only (better to use entire face but we do not have it) ---*/
for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
flux[iVar] = field(iPoint, iVar) / volume;
fluxReflected[iVar] = flux[iVar];
}
/*--- project the flux ---*/
su2double ProjFlux = 0.0;
for (size_t iDim = 0; iDim < nDim; iDim++) ProjFlux += flux[iDim + 1] * UnitNormal[iDim];

for (size_t iDim = 0; iDim < nDim; iDim++)
fluxReflected[iDim + 1] = flux[iDim + 1] - 2.0 * ProjFlux * UnitNormal[iDim];
} // loop over the edges

// /*--- For nodes shared with walls, we can simply add the mirrored contribution. The nonmirrored
// * contributions is added
// * in the routine below. ---*/
// if (nodes->GetSolidBoundary(iPoint)) {
// cout << "point shared by symmetry and wall" << endl;
// su2double volume = nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint);
// const auto area = geometry.vertex[iMarker][iVertex]->GetNormal();
// const su2double* VertexNormal = geometry.vertex[iMarker][iVertex]->GetNormal();
// const auto NormArea = GeometryToolbox::Norm(nDim, VertexNormal);
// su2double UnitNormal[nDim] = {0.0};
// for (size_t iDim = 0; iDim < nDim; iDim++) UnitNormal[iDim] = VertexNormal[iDim] / NormArea;
// su2double ProjArea = 0.0;
// for (unsigned long iDim = 0; iDim < nDim; iDim++) ProjArea += area[iDim] * UnitNormal[iDim];
// su2double areaReflected[nDim] = {0.0};
// for (size_t iDim = 0; iDim < nDim; iDim++)
// areaReflected[iDim] = area[iDim] - 2.0 * ProjArea * UnitNormal[iDim];

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

// /*--- First, use the values at node i only (better to use entire face but we do not have it) ---*/
// for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
// flux[iVar] = field(iPoint, iVar) / volume;
// fluxReflected[iVar] = flux[iVar];
// }
// /*--- project the flux ---*/
// su2double ProjFlux = 0.0;
// for (size_t iDim = 0; iDim < nDim; iDim++) ProjFlux += flux[iDim + 1] * UnitNormal[iDim];

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

// for (size_t iDim = 0; iDim < nDim; iDim++)
// fluxReflected[iDim + 1] = flux[iDim + 1] - 2.0 * ProjFlux * UnitNormal[iDim];

// for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
// for (size_t iDim = 0; iDim < nDim; ++iDim) {
// gradient(iPoint, iVar, iDim) -= 0.5*fluxReflected[iVar] * areaReflected[iDim];
// gradient(iPoint, iVar, iDim) += 0.5*flux[iVar] * area[iDim];
// }
// }
// }

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

// // // /*--- check if point is shared with an inlet, outlet or far_field ---*/
// if (nodes->Getinoutfar(iPoint)) {
// cout << "point = " << iPoint << endl;
// su2double volume = nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint);
// const auto area = geometry.vertex[iMarker][iVertex]->GetNormal();
// const su2double* VertexNormal = geometry.vertex[iMarker][iVertex]->GetNormal();
// const auto NormArea = GeometryToolbox::Norm(nDim, VertexNormal);
// su2double UnitNormal[nDim] = {0.0};
// for (size_t iDim = 0; iDim < nDim; iDim++) UnitNormal[iDim] = VertexNormal[iDim] / NormArea;
// su2double ProjArea = 0.0;
// for (unsigned long iDim = 0; iDim < nDim; iDim++) ProjArea += area[iDim] * UnitNormal[iDim];
// su2double areaReflected[nDim] = {0.0};
// for (size_t iDim = 0; iDim < nDim; iDim++)
// areaReflected[iDim] = area[iDim] - 2.0 * ProjArea * UnitNormal[iDim];

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

// /*--- First, use the values at node i only (better to use entire face but we do not have it) ---*/
// for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
// // note that we have 2x the regular volume
// // we use 1x the regular volume in the other routine so we have
// // (1/2)*(1/Volume) * (F*a + Fr*ar) = (1/2)*(1/Volume)*(Fr*ar) + 1/volume*(F*a) - (1/2)*
// flux[iVar] = field(iPoint, iVar) / volume;
// fluxReflected[iVar] = flux[iVar];
// }
// /*--- project the flux ---*/
// su2double ProjFlux = 0.0;
// for (size_t iDim = 0; iDim < nDim; iDim++) ProjFlux += flux[iDim + 1] * UnitNormal[iDim];

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

// for (size_t iDim = 0; iDim < nDim; iDim++)
// fluxReflected[iDim + 1] = flux[iDim + 1] - 2.0 * ProjFlux * UnitNormal[iDim];

// for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
// for (size_t iDim = 0; iDim < nDim; ++iDim) {
// cout << "before: " << iPoint << ", " << iVar << ", " << iDim << ", " << gradient(iPoint,iVar,iDim)
// << ", delta_L: " << flux[iVar] <<", "<< area[iDim]
// << ", delta_R: " << fluxReflected[iVar] <<", "<< areaReflected[iDim]
// << ", " << gradient(iPoint,iVar,iDim) - fluxReflected[iVar] * areaReflected[iDim] << endl;
// gradient(iPoint, iVar, iDim) -= 0.5*fluxReflected[iVar] * areaReflected[iDim];
// // below we subtract 1/V, but we need only 0.5/V, so we add 0.5*V here
// gradient(iPoint, iVar, iDim) += 0.5*flux[iVar] * area[iDim];
// }
// }
// }

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
for (size_t iDim = 0; iDim < nDim; ++iDim) {
gradient(iPoint, iVar, iDim) -= fluxReflected[iVar] * areaReflected[iDim];
}
}
}

}

} //ivertex
} //symmetry
Expand All @@ -262,18 +313,89 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
if (!nodes->GetDomain(iPoint)) continue;

su2double volume = nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint);

const auto area = geometry.vertex[iMarker][iVertex]->GetNormal();
for (size_t iVar = varBegin; iVar < varEnd; iVar++)
flux[iVar] = field(iPoint,iVar) / volume;

// 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;

// we have to find the edges that were missing in the symmetry computations.
// So we find the jPoints that are on the inlet plane
// so we loop over all neighbor of iPoint, find all jPoints and then check if it is on the inlet
for (size_t iNeigh = 0; iNeigh < nodes->GetnPoint(iPoint); ++iNeigh) {
size_t iEdge = nodes->GetEdge(iPoint, iNeigh);

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable iEdge is not used.
size_t jPoint = nodes->GetPoint(iPoint, iNeigh);
if (nodes->Getinoutfar(jPoint)) {
cout << " jPoint " << jPoint << " is on the inlet plane" << endl;
// this edge jPoint - jPoint is the missing edge for the symmetry computations
//compute the flux on the face between iPoint and jPoint
//for (size_t iVar = varBegin; iVar < varEnd; iVar++) {
// flux[iVar] = 0.5*(field(iPoint,iVar) + field(jPoint, iVar)) / (2.0*volume);
}
if (nodes->GetSymmetry(jPoint)) {
cout << " jPoint " << jPoint << " is on the symmetry plane" << endl;
// this edge iPoint - jPoint is the missing edge for the symmetry computations
// we now need to get the normal of the symmetry plane at jpoint.
// so we loop over the markers, find all symmetry planes, check if the ipoint is on the plane
const su2double* VertexNormal = getVertexNormalfromPoint(config, geometry,jPoint);
cout << " vertex normal = " << VertexNormal[0] <<", " << VertexNormal[1] << endl;
// get the normal on the vertex

// now reflect in the mirror
// reflected normal V=U - 2U_t
const auto NormArea = GeometryToolbox::Norm(nDim, VertexNormal);
su2double UnitNormal[nDim] = {0.0};
for (size_t iDim = 0; iDim < nDim; iDim++)
UnitNormal[iDim] = VertexNormal[iDim] / NormArea;
su2double ProjArea = 0.0;
for (unsigned long iDim = 0; iDim < nDim; iDim++)
ProjArea += area[iDim] * UnitNormal[iDim];
su2double areaReflected[nDim] = {0.0};
for (size_t iDim = 0; iDim < nDim; iDim++)
areaReflected[iDim] = area[iDim] - 2.0 * ProjArea * UnitNormal[iDim];

for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
flux[iVar] = 0.5 * (field(iPoint, iVar) + field(jPoint, iVar)) / (2.0*volume);
fluxReflected[iVar] = flux[iVar];
}

for (size_t iVar = varBegin; iVar < varEnd; iVar++) {
su2double flux = field(iPoint, iVar) / volume;
su2double ProjFlux = 0.0;
for (size_t iDim = 0; iDim < nDim; iDim++)
ProjFlux += flux[iDim + 1] * UnitNormal[iDim];

for (size_t iDim = 0; iDim < nDim; iDim++) gradient(iPoint, iVar, iDim) -= flux * area[iDim];
}
}
for (size_t iDim = 0; iDim < nDim; iDim++)
fluxReflected[iDim + 1] = flux[iDim + 1] - 2.0 * ProjFlux * UnitNormal[iDim];

for (size_t iVar = varBegin; iVar < varEnd; ++iVar) {
for (size_t iDim = 0; iDim < nDim; ++iDim) {
gradient(iPoint, iVar, iDim) += 0.5 * (flux[iVar] * area[iDim] + fluxReflected[iVar] *
areaReflected[iDim]);
}
}

} // if symmetry
} //neighbors

} else {
// if we are on a marker but not on a share point between a symmetry and an inlet/outlet
for (size_t iVar = varBegin; iVar < varEnd; iVar++) {
for (size_t iDim = 0; iDim < nDim; iDim++) {
gradient(iPoint, iVar, iDim) -= flux[iVar] * area[iDim];
}
} // loop over variables
} // symmetry and in/out shared node
} // vertices
END_SU2_OMP_FOR
}
}
} //found right marker
} // iMarkers

/*--- If no solver was provided we do not communicate ---*/

Expand All @@ -293,6 +415,8 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
}
} // namespace detail



/*!
* \brief Instantiations for 2D and 3D.
* \ingroup FvmAlgos
Expand Down

0 comments on commit 7387f8c

Please sign in to comment.