From 7387f8c22aaa1f810ce0c8c53fc94a9f830773aa Mon Sep 17 00:00:00 2001 From: bigfooted Date: Mon, 8 Jan 2024 23:16:26 +0100 Subject: [PATCH] corrected GG gradient computation on shared nodes --- Common/include/CConfig.hpp | 7 + Common/include/geometry/dual_grid/CPoint.hpp | 3 + Common/src/CConfig.cpp | 8 +- Common/src/geometry/CPhysicalGeometry.cpp | 8 +- Common/src/geometry/dual_grid/CPoint.cpp | 1 + .../gradients/computeGradientsGreenGauss.hpp | 244 +++++++++++++----- 6 files changed, 206 insertions(+), 65 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 640514c0067..10a714b9523 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -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 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 diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index 55dc1d6f7a1..2c3d3f38a78 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -75,6 +75,7 @@ class CPoint { su2vector SolidBoundary; /*!< \brief To see if a point belong to the physical boundary (without includin MPI). */ su2vector inoutfar; + su2vector Symmetry; su2vector ViscousBoundary; /*!< \brief To see if a point belong to the physical boundary (without includin MPI). */ su2vector @@ -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. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 01ebe366a86..36810830443 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -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) { diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 2176d175f9b..6aa120fb2c2 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -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 "<GetMarker_All_KindBC(iMarker) << endl; if (config->Getinoutfar(iMarker)) {cout << "inoutfar + // found."<Setinoutfar(Point_Surface, true);} if (config->GetSymmetry(iMarker)) {cout << "symmetry + // found"<SetSymmetry(Point_Surface, true);} if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) nodes->SetPeriodicBoundary(Point_Surface, true); } @@ -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); diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index db81b18a365..2b8a908b177 100644 --- a/Common/src/geometry/dual_grid/CPoint.cpp +++ b/Common/src/geometry/dual_grid/CPoint.cpp @@ -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; diff --git a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp index 90b617d479c..987041db079 100644 --- a/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp +++ b/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp @@ -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 @@ -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 + su2double* jcoord = nodes->GetCoord(jPoint); // nijso: debugging - 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); /*--- Determine if edge points inwards or outwards of iPoint. * If inwards we need to flip the area vector. ---*/ @@ -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]; + + // /*--- 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) -= 0.5*fluxReflected[iVar] * areaReflected[iDim]; + // gradient(iPoint, iVar, iDim) += 0.5*flux[iVar] * area[iDim]; + // } + // } + // } + + // // // /*--- 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]; + + // /*--- 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]; + + // 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]; + // } + // } + // } - 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 @@ -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); + 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 ---*/ @@ -293,6 +415,8 @@ void computeGradientsGreenGauss(CSolver* solver, MPI_QUANTITIES kindMpiComm, PER } } // namespace detail + + /*! * \brief Instantiations for 2D and 3D. * \ingroup FvmAlgos