Skip to content

Commit

Permalink
fig GG share nodes interpolation
Browse files Browse the repository at this point in the history
  • Loading branch information
@bigfooted
bigfooted committed Jan 15, 2024
1 parent 8e346b0 commit e9cb623
Showing 1 changed file with 10 additions and 37 deletions.
47 changes: 10 additions & 37 deletions SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,7 +150,6 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
for (size_t iVertex = 0; iVertex < geometry.GetnVertex(iMarker); ++iVertex) {

size_t iPoint = geometry.vertex[iMarker][iVertex]->GetNode();
//cout << "iPoint = " << iPoint << endl;
auto nodes = geometry.nodes;
// we need to set the gradient to zero for the entire marker to prevent double-counting
// points that are shared by other markers
Expand All @@ -159,13 +158,8 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER

su2double halfOnVol = 0.5 / (nodes->GetVolume(iPoint) + nodes->GetPeriodicVolume(iPoint));
for (size_t iNeigh = 0; iNeigh < nodes->GetnPoint(iPoint); ++iNeigh) {
//cout << " edge = "<<iNeigh << "/"<<nodes->GetnPoint(iPoint) << endl;
size_t iEdge = nodes->GetEdge(iPoint, iNeigh);
//cout <<"edge " << iEdge << endl;
size_t jPoint = nodes->GetPoint(iPoint, iNeigh);
/*su2double* icoord = nodes->GetCoord(iPoint); // nijso: debugging */
/*su2double* jcoord = nodes->GetCoord(jPoint); // nijso: debugging */
//cout << " jPoint="<<jPoint << endl;

/*--- Determine if edge points inwards or outwards of iPoint.
* If inwards we need to flip the area vector. ---*/
Expand Down Expand Up @@ -237,51 +231,36 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER

// 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)) ||
(nodes->GetSymmetry(iPoint) && nodes->GetSolidBoundary(iPoint)) ) {
//cout << "iPoint " << iPoint <<" is on the inlet and on the symmetry" << endl;
//cout << "coords = " << nodes->GetCoord(iPoint)[0] << ", " <<nodes->GetCoord(iPoint)[1] << ", " << nodes->GetCoord(iPoint)[2] << 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
if ( (nodes->GetSymmetry(iPoint) && nodes->Getinoutfar(iPoint)) ||
(nodes->GetSymmetry(iPoint) && nodes->GetSolidBoundary(iPoint)) ) {
// 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
unsigned long jPoint = 0;
for (size_t iNeigh = 0; iNeigh < nodes->GetnPoint(iPoint); ++iNeigh) {
jPoint = nodes->GetPoint(iPoint, iNeigh);
if (nodes->Getinoutfar(jPoint) || nodes->GetSolidBoundary(jPoint)) {
//cout << " jPoint " << jPoint << " is on the inlet plane" << endl;
//if (nodes->GetSymmetry(jPoint))
// cout << "jpoint is on the symmetry" << endl;
//else
// cout << "jpoint is not on the symmetry" << endl;

//cout << "coords = " << nodes->GetCoord(jPoint)[0] << ", " <<nodes->GetCoord(jPoint)[1] << ", " << nodes->GetCoord(jPoint)[2] << endl;


su2double dir = (iPoint < jPoint) ? 1.0 : -1.0;
su2double dir = 1.0; //(iPoint < jPoint) ? 1.0 : -1.0;
su2double weight = dir / (2.0*volume);

// 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) {
// average of i and j at the midway point
flux[iVar] = weight * 0.5 * (field(iPoint, iVar) + field(jPoint, iVar));
/* original : flux[iVar] = field(iPoint,iVar) / volume; */
// average on the face between iPoint and the midway point on the dual edge.
flux[iVar] = weight * (0.75 * field(iPoint, iVar) + 0.25 *field(jPoint, iVar));
fluxReflected[iVar] = flux[iVar];
//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;
// 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;
//cout << " vertex unit normal = " << UnitNormal[0] <<", " << UnitNormal[1] << endl;

su2double ProjArea = 0.0;
for (unsigned long iDim = 0; iDim < nDim; iDim++)
Expand All @@ -291,28 +270,22 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER
for (size_t iDim = 0; iDim < nDim; iDim++)
areaReflected[iDim] = area[iDim] - 2.0 * ProjArea * UnitNormal[iDim];

// velocity components of the flux
su2double ProjFlux = 0.0;
for (size_t iDim = 0; iDim < nDim; iDim++)
ProjFlux += flux[iDim + 1] * UnitNormal[iDim];
// velocity components of the flux
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) {
//cout << "fluxReflected ("<<iVar<<")="<<fluxReflected[iVar]<<endl;

for (size_t iDim = 0; iDim < nDim; ++iDim) {
gradient(iPoint, iVar, iDim) += (flux[iVar] * area[iDim] + fluxReflected[iVar]*areaReflected[iDim]) ;
//cout << "iDim="<<iDim << ", " << -flux[iVar] * area[iDim] << " " << fluxReflected[iVar] * areaReflected[iDim]<< " " << gradient(iPoint, iVar, iDim) << endl;

gradient(iPoint, iVar, iDim) -= (flux[iVar] * area[iDim] + fluxReflected[iVar]*areaReflected[iDim]) ;
}
}
}
}

} else {
} else {
for (size_t iVar = varBegin; iVar < varEnd; iVar++)
flux[iVar] = field(iPoint,iVar) / volume;

Expand Down

0 comments on commit e9cb623

Please sign in to comment.