diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 9d3e574b4328..1adcd124e4b7 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -62,6 +62,7 @@ extern "C" { #include "../CConfig.hpp" #include "../toolboxes/graph_toolbox.hpp" +#include "../toolboxes/geometry_toolbox.hpp" #include "../adt/CADTElemClass.hpp" using namespace std; @@ -517,6 +518,32 @@ class CGeometry { */ void AllocatePeriodicComms(unsigned short val_countPerPeriodicPoint); + /*! \brief Fraction of an edge stencil owned by this periodic copy. */ + inline passivedouble GetPeriodicEdgeWeight(unsigned long iPoint, unsigned long jPoint, const CConfig& config) const { + if (!nodes->GetPeriodicBoundary(iPoint) || !nodes->GetPeriodicBoundary(jPoint)) return 1.0; + passivedouble weight = 1.0; + for (auto iMarker = 0u; iMarker < nMarker; ++iMarker) { + if (config.GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY || nodes->GetVertex(iPoint, iMarker) < 0 || + nodes->GetVertex(jPoint, iMarker) < 0) + continue; + const auto tag = config.GetMarker_All_TagBound(iMarker); + const auto donor = config.GetMarker_Periodic_Donor(tag); + if (donor < nMarker && nodes->GetVertex(iPoint, donor) >= 0 && nodes->GetVertex(jPoint, donor) >= 0) { + /*--- An edge along the rotation axis is shared by all N sectors, + * rather than by two copies of each of the two periodic faces. ---*/ + if (config.GetMarker_All_PerBound(iMarker) > config.GetnMarker_Periodic() / 2) continue; + const auto* angles = config.GetPeriodicRotAngles(tag); + su2double rotation[3][3]; + GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation); + const auto cosAngle = SU2_TYPE::GetValue(0.5 * (rotation[0][0] + rotation[1][1] + rotation[2][2] - 1)); + weight *= acos(std::max(-1.0, std::min(1.0, cosAngle))) / (2 * PI_NUMBER); + } else { + weight *= 0.5; + } + } + return weight; + } + /*! * \brief Routine to launch non-blocking recvs only for all periodic communication with neighboring partitions. * \note This routine is called by any class that has loaded data into the generic communication buffers. diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index d3bca729bf5b..e02d613962ec 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1157,7 +1157,14 @@ void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* co auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); if (fine_grid->nodes->GetBoundary(iFinePoint)) { nodes->SetBoundary(iCoarsePoint, nMarker); - break; + nodes->SetPeriodicBoundary(iCoarsePoint, nodes->GetPeriodicBoundary(iCoarsePoint) || + fine_grid->nodes->GetPeriodicBoundary(iFinePoint)); + nodes->SetPhysicalBoundary(iCoarsePoint, nodes->GetPhysicalBoundary(iCoarsePoint) || + fine_grid->nodes->GetPhysicalBoundary(iFinePoint)); + nodes->SetSolidBoundary( + iCoarsePoint, nodes->GetSolidBoundary(iCoarsePoint) || fine_grid->nodes->GetSolidBoundary(iFinePoint)); + nodes->SetViscousBoundary( + iCoarsePoint, nodes->GetViscousBoundary(iCoarsePoint) || fine_grid->nodes->GetViscousBoundary(iFinePoint)); } } diff --git a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp index 296f433ff40f..a239f62989d6 100644 --- a/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp +++ b/SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp @@ -253,6 +253,8 @@ void computeGradientsLeastSquares(CSolver* solver, if (weight > 0.0) { weight = 1.0 / weight; + if (config.GetnMarker_Periodic() > 2) + weight *= geometry.GetPeriodicEdgeWeight(iPoint, jPoint, config); for (size_t iDim = 0; iDim < nDim; ++iDim) for (size_t jDim = iDim; jDim < nDim; ++jDim) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index ae07bef84552..10e0bf538cbd 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -807,8 +807,10 @@ class CFVMFlowSolverBase : public CSolver { iPoint_UndLapl[iPoint] = fmax(iPoint_UndLapl[iPoint], fabs(sensVar_j - sensVar_i) / fmin(sensVar_j, sensVar_i)); } else { /*--- Jameson dissipation sensor, add variable difference and variable sum. ---*/ - iPoint_UndLapl[iPoint] += sensVar_j - sensVar_i; - jPoint_UndLapl[iPoint] += sensVar_j + sensVar_i; + const auto weight = config->GetnMarker_Periodic() > 2 ? + geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) : 1.0; + iPoint_UndLapl[iPoint] += weight * (sensVar_j - sensVar_i); + jPoint_UndLapl[iPoint] += weight * (sensVar_j + sensVar_i); } } diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index cb44757c6280..b68c9de9e5d2 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -260,6 +260,13 @@ void CFVMFlowSolverBase::CommunicateInitialState(CGeometry* geometry, cons InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); } + + /*--- The periodic communication updates the number of neighbors of the owned points only, update the halos. ---*/ + + if (config->GetnMarker_Periodic() > 0) { + geometry->InitiateComms(geometry, config, MPI_QUANTITIES::NEIGHBORS); + geometry->CompleteComms(geometry, config, MPI_QUANTITIES::NEIGHBORS); + } SetImplicitPeriodic(euler_implicit); if (MGLevel == MESH_0) SetRotatePeriodic(true); diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 7a7b6c726021..feda8a8d6bb2 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -133,6 +133,8 @@ class CSolver { /*--- End variables that need to go. ---*/ + su2activevector periodicNeighborCount; /*!< \brief Partial periodic stencil count, including shared edges. */ + su2activevector iPoint_UndLapl; /*!< \brief Auxiliary variable for the undivided Laplacians. */ su2activevector jPoint_UndLapl; /*!< \brief Auxiliary variable for the undivided Laplacians. */ diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 32c265e55656..5f056d860473 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2404,13 +2404,16 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *co /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ if (boundary_i && !boundary_j) continue; + const auto weight = config->GetnMarker_Periodic() > 2 ? + geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) : 1.0; + /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ for (unsigned short iVar = 0; iVar < nVar; iVar++) - nodes->AddUnd_Lapl(iPoint, iVar, nodes->GetSolution(jPoint,iVar)-nodes->GetSolution(iPoint,iVar)); + nodes->AddUnd_Lapl(iPoint, iVar, weight * (nodes->GetSolution(jPoint,iVar)-nodes->GetSolution(iPoint,iVar))); su2double Pressure_j = nodes->GetPressure(jPoint); - nodes->AddUnd_Lapl(iPoint, nVar-1, Pressure_j-Pressure_i); + nodes->AddUnd_Lapl(iPoint, nVar-1, weight * (Pressure_j-Pressure_i)); } } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index e23d25af07c2..2aacbf2cc4b3 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -148,12 +148,16 @@ CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, const CSolver* fl ghostNodes = make_unique(Solution_Inf[0], maxMarkerVertices, nDim, nVar, config); } - /*--- Communicate and store volume and the number of neighbors for any dual CVs that lie on on periodic markers. ---*/ - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); + /*--- Communicate and store volume and the number of neighbors for any dual CVs that lie on on periodic markers. + * With a flow solver on the same geometry this was already done by the flow solver, and the values are + * accumulated, so it must not be done twice. ---*/ + if (!flow) { + for (auto iPeriodic = 1u; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_VOLUME); + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_NEIGHBORS); + } } /*--- Store if implicit scheme is used. This has implications on the Residual and Jacobian handling for periodic * boundaries ---*/ diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 56cee4b6ad52..895da0e329e4 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -207,7 +207,7 @@ void CSolver::GetPeriodicCommCountAndType(const CConfig* config, break; case PERIODIC_NEIGHBORS: COUNT_PER_POINT = 1; - MPI_TYPE = COMM_TYPE::UNSIGNED_SHORT; + MPI_TYPE = COMM_TYPE::DOUBLE; break; case PERIODIC_RESIDUAL: COUNT_PER_POINT = nVar + nVar*nVar + 1; @@ -352,6 +352,18 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, SU2_MPI::Error("The NEMO solvers do not support rotational periodicity yet.", CURRENT_FUNCTION); } + const bool intersectingPairs = config->GetnMarker_Periodic() > 2; + if (commType == PERIODIC_NEIGHBORS && intersectingPairs && val_periodic_index == 1) { + SU2_OMP_SAFE_GLOBAL_ACCESS(periodicNeighborCount.resize(geometry->GetnPoint());) + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) { + periodicNeighborCount(iPoint) = 0; + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) + periodicNeighborCount(iPoint) += geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config); + } + END_SU2_OMP_FOR + } + /*--- Local variables ---*/ bool boundary_i, boundary_j; @@ -413,7 +425,6 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, su2double *bufDSend = geometry->bufD_PeriodicSend; - unsigned short *bufSSend = geometry->bufS_PeriodicSend; /*--- Handle the different types of gradient and limiter. ---*/ @@ -501,6 +512,11 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, case PERIODIC_NEIGHBORS: + if (intersectingPairs) { + bufDSend[buf_offset] = periodicNeighborCount(iPoint); + break; + } + nNeighbor = 0; for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { @@ -508,13 +524,13 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, that we avoid double counting neighbors on both sides. If not, increment the count of neighbors for the donor. ---*/ - if (!geometry->nodes->GetPeriodicBoundary(jPoint)) + if (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) == 1.0) nNeighbor++; } /*--- Store the number of neighbors in bufffer. ---*/ - bufSSend[buf_offset] = nNeighbor; + bufDSend[buf_offset] = nNeighbor; break; @@ -605,6 +621,16 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, case PERIODIC_LAPLACIAN: + if (intersectingPairs) { + for (auto iField = 0u; iField < nVar; ++iField) + bufDSend[buf_offset + iField] = base_nodes->GetUndivided_Laplacian(iPoint, iField); + if (rotate_periodic) Rotate(zeros, &bufDSend[buf_offset + 1], &Und_Lapl[1]); + if (rotate_periodic) + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) + bufDSend[buf_offset + 1 + iCoordinate] = Und_Lapl[1 + iCoordinate]; + break; + } + /*--- For JST, the undivided Laplacian must be computed consistently by using the complete control volume info from both sides of the periodic face. ---*/ @@ -617,7 +643,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Avoid periodic boundary points so that we do not duplicate edges on both sides of the periodic BC. ---*/ - if (!geometry->nodes->GetPeriodicBoundary(jPoint)) { + if (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) == 1.0) { /*--- Solution differences ---*/ @@ -671,6 +697,11 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, break; case PERIODIC_SENSOR: { + if (intersectingPairs) { + bufDSend[buf_offset] = iPoint_UndLapl(iPoint); + bufDSend[buf_offset + 1] = jPoint_UndLapl(iPoint); + break; + } const bool msw = config->GetKind_Upwind_Flow() == UPWIND::MSW; /*--- For the centered schemes, the sensor must be computed @@ -683,7 +714,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Avoid halos and boundary points so that we don't duplicate edges on both sides of the periodic BC. ---*/ - if (geometry->nodes->GetPeriodicBoundary(jPoint)) continue; + if (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) < 1.0) continue; /*--- Use density instead of pressure for incomp. flows. ---*/ @@ -703,7 +734,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, if ((!boundary_i || boundary_j) && geometry->nodes->GetDomain(iPoint)) { if (msw) { - Sensor_i = fmax(Sensor_i, fabs(Pressure_j - Pressure_i)) / fmin(Pressure_i, Pressure_j); + Sensor_i = fmax(Sensor_i, fabs(Pressure_j - Pressure_i) / fmin(Pressure_i, Pressure_j)); } else { Sensor_i += (Pressure_j - Pressure_i); Sensor_j += (Pressure_i + Pressure_j); @@ -781,6 +812,42 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Set a flag for unweighted or weighted least-squares. ---*/ + if (intersectingPairs) { + /*--- Rmatrix stores the upper triangle of the normal equations; + * entry (2,1) is a duplicate of (0,2) used by the LS factorization. ---*/ + su2double matrix[3][3] = {}, rotated[3][3] = {}; + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) + for (auto jDim = 0u; jDim < nDim; ++jDim) + matrix[iCoordinate][jDim] = + base_nodes->GetRmatrix(iPoint, min(iCoordinate,jDim), max(iCoordinate,jDim)); + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) + for (auto jDim = 0u; jDim < nDim; ++jDim) + for (auto kDim = 0u; kDim < nDim; ++kDim) + for (auto lDim = 0u; lDim < nDim; ++lDim) { + const auto qik = nDim == 2 ? rotMatrix2D[iCoordinate][kDim] : rotMatrix3D[iCoordinate][kDim]; + const auto qjl = nDim == 2 ? rotMatrix2D[jDim][lDim] : rotMatrix3D[jDim][lDim]; + rotated[iCoordinate][jDim] += qik * matrix[kDim][lDim] * qjl; + } + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) + for (auto jDim = 0u; jDim < nDim; ++jDim) + bufDSend[buf_offset++] = iCoordinate <= jDim ? rotated[iCoordinate][jDim] : + (nDim == 3 && iCoordinate == 2 && jDim == 1 ? rotated[0][2] : su2double(0)); + for (auto iField = 0u; iField < ICOUNT; ++iField) + Rotate(zeros, gradient[iPoint][iField], rotBlock[iField]); + if (rotate_periodic) { + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) { + su2double velocity[3] = {}, rotatedVelocity[3] = {}; + for (auto jDim = 0u; jDim < nDim; ++jDim) velocity[jDim] = rotBlock[1+jDim][iCoordinate]; + Rotate(zeros, velocity, rotatedVelocity); + for (auto jDim = 0u; jDim < nDim; ++jDim) rotBlock[1+jDim][iCoordinate] = rotatedVelocity[jDim]; + } + } + for (auto iField = 0u; iField < ICOUNT; ++iField) + for (auto iCoordinate = 0u; iCoordinate < nDim; ++iCoordinate) + bufDSend[buf_offset++] = rotBlock[iField][iCoordinate]; + break; + } + switch(commType) { case PERIODIC_SOL_ULS: case PERIODIC_SOL_ULS_R: @@ -826,7 +893,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Avoid periodic boundary points so that we do not duplicate edges on both sides of the periodic BC. ---*/ - if (!geometry->nodes->GetPeriodicBoundary(jPoint)) { + if (geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config) == 1.0) { /*--- Get coordinates for the neighbor point. ---*/ @@ -1041,7 +1108,6 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, const su2double *bufDRecv = geometry->bufD_PeriodicRecv; - const unsigned short *bufSRecv = geometry->bufS_PeriodicRecv; /*--- Handle the different types of gradient and limiter. ---*/ @@ -1118,13 +1184,13 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, break; case PERIODIC_NEIGHBORS: - - /*--- Store the extra neighbors on the periodic face. ---*/ - - nNeighbor = (geometry->nodes->GetnNeighbor(iPoint) + - bufSRecv[buf_offset]); - geometry->nodes->SetnNeighbor(iPoint, nNeighbor); - + if (config->GetnMarker_Periodic() > 2) { + periodicNeighborCount(iPoint) += bufDRecv[buf_offset]; + geometry->nodes->SetnNeighbor(iPoint, SU2_TYPE::Int(periodicNeighborCount(iPoint) + 0.5)); + } else { + nNeighbor = geometry->nodes->GetnNeighbor(iPoint) + SU2_TYPE::Int(bufDRecv[buf_offset]); + geometry->nodes->SetnNeighbor(iPoint, nNeighbor); + } break; case PERIODIC_RESIDUAL: @@ -2280,6 +2346,7 @@ void CSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config) for (unsigned short iVar = 0; iVar < nVar; iVar++) { su2double delta = base_nodes->GetSolution(jPoint,iVar)-base_nodes->GetSolution(iPoint,iVar); + if (config->GetnMarker_Periodic() > 2) delta *= geometry->GetPeriodicEdgeWeight(iPoint, jPoint, *config); base_nodes->AddUnd_Lapl(iPoint, iVar, delta); } } diff --git a/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg new file mode 100644 index 000000000000..6ab487b08981 --- /dev/null +++ b/TestCases/incomp_navierstokes/streamwise_periodic/chtPinArray_2d/periodic_weak_heat.cfg @@ -0,0 +1,79 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: 2D pin array, periodic flow driven by a body force, with % +% the weakly coupled heat equation (flow and heat solver on % +% the same mesh, the heat solver does not change the flow) % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +SOLVER= INC_NAVIER_STOKES +KIND_TURB_MODEL= NONE +RESTART_SOL= NO +% +INC_ENERGY_EQUATION= NO +WEAKLY_COUPLED_HEAT_EQUATION= YES + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= CONSTANT +INC_DENSITY_INIT= 1045.0 +INC_VELOCITY_INIT= ( 0.0001, 0.0, 0.0 ) +INC_TEMPERATURE_INIT= 338.0 +INC_NONDIM= DIMENSIONAL +% +SPECIFIC_HEAT_CP= 3540.0 +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 0.001385 +CONDUCTIVITY_MODEL= CONSTANT_PRANDTL +PRANDTL_LAM= 11.7 +% +BODY_FORCE= YES +BODY_FORCE_VECTOR= ( 0.5, 0.0, 0.0 ) + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_SYM= ( fluid_symmetry ) +MARKER_PERIODIC= ( fluid_inlet, fluid_outlet, 0.0,0.0,0.0, 0.0,0.0,0.0, 0.0111544,0.0,0.0 ) +MARKER_HEATFLUX= ( fluid_pin1_interface, 0.0, fluid_pin2_interface, 0.0, fluid_pin3_interface, 0.0 ) +% +MARKER_MONITORING= ( fluid_pin2_interface ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 1000 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1e-8 +LINEAR_SOLVER_ITER= 20 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% -------------------- HEAT NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_HEAT= SCALAR_UPWIND +TIME_DISCRE_HEAT= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= RMS_PRESSURE +CONV_RESIDUAL_MINVAL= -16 +CONV_STARTITER= 10 +ITER= 2000 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= fluid.su2 +OUTPUT_WRT_FREQ= 9999 +SCREEN_OUTPUT= ( INNER_ITER, RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, DRAG ) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index e79219efbeba..998dcc44eeb3 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -788,6 +788,14 @@ def main(): inc_heatTransfer_BC.test_vals = [-8.904266, -7.745636, -8.064003, -1.079197, -1671.000000] test_list.append(inc_heatTransfer_BC) + # 2D pin array, periodic with a body force, flow and weakly coupled heat equation + inc_periodic_weak_heat = TestCase('inc_periodic_weak_heat') + inc_periodic_weak_heat.cfg_dir = "incomp_navierstokes/streamwise_periodic/chtPinArray_2d" + inc_periodic_weak_heat.cfg_file = "periodic_weak_heat.cfg" + inc_periodic_weak_heat.test_iter = 10 + inc_periodic_weak_heat.test_vals = [-4.712559, -5.747243, -5.751727, 873.702714] + test_list.append(inc_periodic_weak_heat) + ############################ ### Incompressible RANS ### ############################ @@ -1172,7 +1180,7 @@ def main(): sbs_backward_step.cfg_dir = "backscatter/backward_step" sbs_backward_step.cfg_file = "backwardStep.cfg" sbs_backward_step.test_iter = 3 - sbs_backward_step.test_vals = [-6.352884, -3.465372, -5.507901, -3.906545, -9.506305, -6.365236, -6.331021, -6.331028] + sbs_backward_step.test_vals = [-6.352885, -3.465372, -5.507901, -3.906544, -9.506305, -6.365236, -6.331021, -6.331028] sbs_backward_step.unsteady = True sbs_backward_step.decompress = True sbs_backward_step.grid_file = "backward_step.su2" diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index 72baf525ca68..671a8e1e262f 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -26,6 +26,8 @@ */ #include "catch.hpp" +#include "../UnitQuadTestCase.hpp" +#include "../../Common/include/geometry/CMultiGridGeometry.hpp" #include "../../Common/include/geometry/CPhysicalGeometry.hpp" #include "../../Common/include/containers/container_decorators.hpp" #include "../../SU2_CFD/include/solvers/CSolver.hpp" @@ -154,3 +156,194 @@ TEST_CASE("GG", "[Gradients]") { testGreenGauss(); } TEST_CASE("LS", "[Gradients]") { testLeastSquares(false); } TEST_CASE("WLS", "[Gradients]") { testLeastSquares(true); } + +TEST_CASE("Intersecting periodic stencils", "[Periodic][Gradients]") { + const auto nPairs = GENERATE(2u, 3u); + const bool diagonal = GENERATE(false, true); + const auto scheme = GENERATE(0u, 1u, 2u); + const bool msw = scheme == 1, incompressible = scheme == 2; + UnitQuadTestCase field; + const auto start = field.config_options.find("MARKER_HEATFLUX="); + const auto end = field.config_options.find("VISCOSITY_MODEL="); + field.config_options.replace(start, end - start, nPairs == 2 ? "MARKER_CUSTOM= (z_plus,z_minus)\n" : ""); + std::string periodic = "MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0, y_minus,y_plus, 0,0,0, 0,0,0, 0,1,0"; + if (nPairs == 3) periodic += ", z_minus,z_plus, 0,0,0, 0,0,0, 0,0,1"; + field.AddOption(periodic + ")"); + field.AddOption(msw ? "CONV_NUM_METHOD_FLOW= MSW" : "CONV_NUM_METHOD_FLOW= JST"); + field.AddOption("MUSCL_FLOW= NO\nNUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES"); + if (incompressible) { + field.SetOption("SOLVER= INC_NAVIER_STOKES"); + field.AddOption("FLUID_MODEL= CONSTANT_DENSITY"); + field.SetOption("KIND_VERIFICATION_SOLUTION= NO_VERIFICATION_SOLUTION"); + } + field.InitConfig(); + field.InitGeometry(true); + if (diagonal) { + /*--- A triangulated plane adds NE/SW edges. The SW donor at a periodic + * corner is reached through two successive pairs, rather than one. ---*/ + std::vector> neighbors(field.geometry->GetnPoint()); + for (auto i = 0ul; i < neighbors.size(); ++i) { + for (auto j : field.geometry->nodes->GetPoints(i)) neighbors[i].push_back(j); + const auto* x = field.geometry->nodes->GetCoord(i); + for (auto j = 0ul; j < neighbors.size(); ++j) { + const auto* y = field.geometry->nodes->GetCoord(j); + if (fabs(y[0] - x[0] - 0.25) < 1e-12 && fabs(y[1] - x[1] - 0.25) < 1e-12 && fabs(y[2] - x[2]) < 1e-12) { + neighbors[i].push_back(j); + neighbors[j].push_back(i); + } + } + } + field.geometry->nodes->SetPoints(neighbors); + delete field.geometry->edges; + field.geometry->SetEdges(); + for (auto i = 0ul; i < neighbors.size(); ++i) field.geometry->nodes->SetnNeighbor(i, neighbors[i].size()); + } + for (auto pair = 1u; pair <= nPairs; ++pair) field.geometry->MatchPeriodic(field.config.get(), pair); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + field.InitSolver(); + auto* solver = field.solver[FLOW_SOL]; + auto* nodes = solver->GetNodes(); + auto value = [](su2double x, su2double y, su2double z) { + return 10 + x * (1 - x) + 2 * y * (1 - y) + 3 * z * (1 - z); + }; + unsigned long target = field.geometry->GetnPoint(); + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + nodes->SetPressure(i, value(x[0], x[1], x[2])); + for (auto v = 0u; v < solver->GetnVar(); ++v) + nodes->SetSolution( + i, v, + v == 0 ? value(x[0], x[1], x[2]) + : (v + 1 == solver->GetnVar() ? (incompressible ? 300.0 : value(x[0], x[1], x[2]) / 0.4) : 0)); + if (field.geometry->nodes->GetDomain(i) && x[0] == 0 && x[1] == 0 && x[2] == (nPairs == 3 ? 0 : 0.5)) target = i; + } + const auto z = nPairs == 3 ? 0.0 : 0.5; + const auto center = value(0, 0, z); + std::vector stencil = {value(0.25, 0, z), value(0.75, 0, z), + value(0, 0.25, z), value(0, 0.75, z), + value(0, 0, z + 0.25), value(0, 0, nPairs == 3 ? 0.75 : z - 0.25)}; + if (diagonal) { + stencil.push_back(value(0.25, 0.25, z)); + stencil.push_back(value(0.75, 0.75, z)); + } + su2double numerator = 0, denominator = 0, maximum = 0; + for (const auto neighbor : stencil) { + numerator += neighbor - center; + denominator += neighbor + center; + maximum = std::max(maximum, fabs(neighbor - center) / std::min(center, neighbor)); + } + field.config->SetGlobalParam(field.config->GetKind_Solver(), RUNTIME_FLOW_SYS); + { + SU2_OMP_PARALLEL { + solver->Preprocessing(field.geometry.get(), field.solver, field.config.get(), 0, 0, RUNTIME_FLOW_SYS, false); + } + END_SU2_OMP_PARALLEL + } + unsigned long localTarget = target < field.geometry->GetnPoint(), globalTarget = 0; + SU2_MPI::Allreduce(&localTarget, &globalTarget, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(globalTarget == 1); + INFO("pairs " << nPairs << ", diagonal " << diagonal << ", MSW " << msw << ", incompressible " << incompressible); + if (target < field.geometry->GetnPoint()) { + CHECK(field.geometry->nodes->GetnNeighbor(target) == stencil.size()); + CHECK(nodes->GetSensor(target) == + Approx(msw ? maximum : (incompressible ? 0.0 : fabs(numerator) / denominator)).margin(1e-12)); + if (!msw) CHECK(nodes->GetUndivided_Laplacian(target, 0) == Approx(numerator).margin(1e-12)); + } + /*--- Periodic scalar field: compare the assembled primitive LS normal + * equations with the complete stencil, including the composed diagonal. ---*/ + solver->SetRotatePeriodic(false); + auto& gradient = nodes->GetGradient_Primitive(); + auto& matrix = nodes->GetRmatrix(); + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v) nodes->SetPrimitive(i, v, value(x[0], x[1], x[2])); + } + SU2_OMP_PARALLEL { + computeGradientsLeastSquares(solver, MPI_QUANTITIES::PRIMITIVE_GRADIENT, PERIODIC_PRIM_LS, *field.geometry, + *field.config, true, nodes->GetPrimitive(), 0, solver->GetnPrimVarGrad(), -1, gradient, + matrix); + } + END_SU2_OMP_PARALLEL + if (target < field.geometry->GetnPoint()) { + CHECK(gradient(target, 0, 0) == Approx(0).margin(1e-12)); + CHECK(gradient(target, 0, 1) == Approx(0).margin(1e-12)); + CHECK(gradient(target, 0, 2) == Approx(nPairs == 3 ? 0 : 3 * (1 - 2 * z)).margin(1e-12)); + } +} + +TEST_CASE("Coarse periodic boundary flags", "[Periodic]") { + UnitQuadTestCase field; + const auto start = field.config_options.find("MARKER_HEATFLUX="); + const auto end = field.config_options.find("VISCOSITY_MODEL="); + field.config_options.replace(start, end - start, + "MARKER_HEATFLUX= (y_minus,0,y_plus,0)\nMARKER_CUSTOM= (z_plus,z_minus)\n"); + field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)"); + field.AddOption("MGLEVEL= 1"); + field.InitConfig(); + field.InitGeometry(true); + CMultiGridGeometry coarse(field.geometry.get(), field.config.get(), 1); + coarse.SetPoint_Connectivity(field.geometry.get()); + coarse.SetVertex(field.geometry.get(), field.config.get()); + unsigned long periodic = 0, missing = 0; + for (auto marker = 0u; marker < coarse.GetnMarker(); ++marker) { + if (field.config->GetMarker_All_KindBC(marker) == HEAT_FLUX) { + for (auto v = 0ul; v < coarse.GetnVertex(marker); ++v) { + const auto i = coarse.vertex[marker][v]->GetNode(); + CHECK(coarse.nodes->GetPhysicalBoundary(i)); + CHECK(coarse.nodes->GetSolidBoundary(i)); + CHECK(coarse.nodes->GetViscousBoundary(i)); + } + } + if (field.config->GetMarker_All_KindBC(marker) != PERIODIC_BOUNDARY) continue; + for (auto vertex = 0ul; vertex < coarse.GetnVertex(marker); ++vertex) { + ++periodic; + if (!coarse.nodes->GetPeriodicBoundary(coarse.vertex[marker][vertex]->GetNode())) ++missing; + } + } + unsigned long local[] = {periodic, missing}, global[2] = {}; + SU2_MPI::Allreduce(local, global, 2, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(global[0] > 0); + CHECK(global[1] == 0); +} + +TEST_CASE("One-cell periodic stencil", "[Periodic][Gradients]") { + UnitQuadTestCase field; + const auto start = field.config_options.find("MARKER_HEATFLUX="); + const auto end = field.config_options.find("VISCOSITY_MODEL="); + field.config_options.replace(start, end - start, "MARKER_CUSTOM= (y_minus,y_plus,z_plus,z_minus)\n"); + field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)"); + field.SetOption("MESH_BOX_SIZE= 2,5,5"); + field.AddOption("NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES"); + field.InitConfig(); + field.InitGeometry(true); + field.geometry->MatchPeriodic(field.config.get(), 1); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + field.InitSolver(); + auto* solver = field.solver[FLOW_SOL]; + auto* nodes = solver->GetNodes(); + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v) nodes->SetPrimitive(i, v, 3 + x[1] + x[2]); + } + solver->SetRotatePeriodic(false); + auto& gradient = nodes->GetGradient_Primitive(); + auto& matrix = nodes->GetRmatrix(); + SU2_OMP_PARALLEL { + computeGradientsLeastSquares(solver, MPI_QUANTITIES::PRIMITIVE_GRADIENT, PERIODIC_PRIM_LS, *field.geometry, + *field.config, true, nodes->GetPrimitive(), 0, solver->GetnPrimVarGrad(), -1, gradient, + matrix); + } + END_SU2_OMP_PARALLEL + unsigned long local = 0, global = 0; + for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + if (x[1] != 0.5 || x[2] != 0.5) continue; + ++local; + CHECK(field.geometry->nodes->GetnNeighbor(i) == 6); + CHECK(matrix(i, 0, 0) == Approx(2).margin(1e-12)); + CHECK(gradient(i, 0, 1) == Approx(1).margin(1e-12)); + CHECK(gradient(i, 0, 2) == Approx(1).margin(1e-12)); + } + SU2_MPI::Allreduce(&local, &global, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(global == 2); +} diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp index d8da75dc0960..d943c1696ee5 100644 --- a/UnitTests/UnitQuadTestCase.hpp +++ b/UnitTests/UnitQuadTestCase.hpp @@ -61,6 +61,16 @@ struct UnitQuadTestCase { */ void AddOption(const std::string& optionLine) { config_options += optionLine + "\n"; } + /*! \brief Replace one existing base option without repeating its key. */ + void SetOption(const std::string& optionLine) { + const auto key = optionLine.substr(0, optionLine.find('=') + 1); + const auto start = config_options.find(key); + if (start == std::string::npos) + AddOption(optionLine); + else + config_options.replace(start, config_options.find('\n', start) - start, optionLine); + } + /*! * \brief Initialize the config structure */ @@ -83,10 +93,13 @@ struct UnitQuadTestCase { /*! * \brief Initialize the geometry */ - void InitGeometry() { + void InitGeometry(bool partition = false) { cout.rdbuf(nullptr); { - auto aux_geometry = std::unique_ptr(new CPhysicalGeometry(config.get(), 0, 1)); + const auto rank = partition ? SU2_MPI::GetRank() : 0; + const auto size = partition ? SU2_MPI::GetSize() : 1; + auto aux_geometry = std::unique_ptr(new CPhysicalGeometry(config.get(), rank, size)); + if (partition) aux_geometry->SetColorGrid_Parallel(config.get()); geometry = std::unique_ptr(new CPhysicalGeometry(aux_geometry.get(), config.get())); } geometry->SetSendReceive(config.get());