From 64d8ff6b50d90d98eb0b2ec979634a4661d595af Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 28 Dec 2024 12:13:32 -0800 Subject: [PATCH] remove axis-aligned restriction --- Common/include/linear_algebra/CSysMatrix.hpp | 4 +- Common/src/linear_algebra/CSysMatrix.cpp | 58 ++++++++++++++------ SU2_CFD/src/solvers/CFEASolver.cpp | 53 ++++++++++++------ TestCases/parallel_regression.py | 2 +- 4 files changed, 78 insertions(+), 39 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 815d52e0708..422a2068a88 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -805,10 +805,10 @@ class CSysMatrix { void EnforceSolutionAtNode(unsigned long node_i, const OtherType* x_i, CSysVector& b); /*! - * \brief Version of EnforceSolutionAtNode for a single degree of freedom. + * \brief Similar to EnforceSolutionAtNode, but for 0 projection in a given direction. */ template - void EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar, OtherType x_i, CSysVector& b); + void EnforceZeroProjection(unsigned long node_i, const OtherType* n, CSysVector& b); /*! * \brief Sets the diagonal entries of the matrix as the sum of the blocks in the corresponding column. diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 4196f638678..9ef0c7b88de 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -1027,36 +1027,59 @@ void CSysMatrix::EnforceSolutionAtNode(const unsigned long node_i, c template template -void CSysMatrix::EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar, OtherType x_i, - CSysVector& b) { +void CSysMatrix::EnforceZeroProjection(unsigned long node_i, const OtherType* n, CSysVector& b) { for (auto index = row_ptr[node_i]; index < row_ptr[node_i + 1]; ++index) { const auto node_j = col_ind[index]; - /*--- Delete row iVar of block j on row i (bij) and ATTEMPT - * to delete column iVar block i on row j (bji). ---*/ + /*--- Remove product components of block j on row i (bij) and ATTEMPT + * to remove solution components of block i on row j (bji). + * This is identical to symmetry correction applied to gradients + * but extended to the entire matrix. ---*/ auto bij = &matrix[index * nVar * nVar]; auto bji = GetBlock(node_j, node_i); - /*--- The "attempt" part. ---*/ + /*--- Attempt to remove solution components. ---*/ + ScalarType nbn{}; if (bji != nullptr) { - for (auto jVar = 0ul; jVar < nVar; ++jVar) { - /*--- Column product. ---*/ - b[node_j * nVar + jVar] -= bji[jVar * nVar + iVar] * x_i; - /*--- Delete entries. ---*/ - bji[jVar * nVar + iVar] = 0.0; + for (auto iVar = 0ul; iVar < nVar; ++iVar) { + ScalarType proj{}; + for (auto jVar = 0ul; jVar < nVar; ++jVar) { + proj += bji[iVar * nVar + jVar] * PassiveAssign(n[jVar]); + } + for (auto jVar = 0ul; jVar < nVar; ++jVar) { + bji[iVar * nVar + jVar] -= proj * PassiveAssign(n[jVar]); + } + nbn += proj * PassiveAssign(n[iVar]); } } - /*--- Delete row. ---*/ - for (auto jVar = 0ul; jVar < nVar; ++jVar) bij[iVar * nVar + jVar] = 0.0; + /*--- Product components. ---*/ + for (auto jVar = 0ul; jVar < nVar; ++jVar) { + ScalarType proj{}; + for (auto iVar = 0ul; iVar < nVar; ++iVar) { + proj += bij[iVar * nVar + jVar] * PassiveAssign(n[iVar]); + } + for (auto iVar = 0ul; iVar < nVar; ++iVar) { + bij[iVar * nVar + jVar] -= proj * PassiveAssign(n[iVar]); + } + } - /*--- Set the diagonal entry of the block to 1. ---*/ - if (node_j == node_i) bij[iVar * (nVar + 1)] = 1.0; + /*--- This part doesn't have the "*2" factor because the product components + * were removed from the result of removing the solution components + * instead of from the original block (bji == bij). ---*/ + if (node_i == node_j) { + for (auto iVar = 0ul; iVar < nVar; ++iVar) { + for (auto jVar = 0ul; jVar < nVar; ++jVar) { + bij[iVar * nVar + jVar] += PassiveAssign(n[iVar]) * nbn * PassiveAssign(n[jVar]); + } + } + } } - /*--- Set known solution in rhs vector. ---*/ - b(node_i, iVar) = x_i; + OtherType proj{}; + for (auto iVar = 0ul; iVar < nVar; ++iVar) proj += b(node_i, iVar) * n[iVar]; + for (auto iVar = 0ul; iVar < nVar; ++iVar) b(node_i, iVar) -= proj * n[iVar]; } template @@ -1203,8 +1226,7 @@ void CSysMatrix::ComputePastixPreconditioner(const CSysVector; \ template void CSysMatrix::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector&); \ - template void CSysMatrix::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double, \ - CSysVector&); \ + template void CSysMatrix::EnforceZeroProjection(unsigned long, const su2double*, CSysVector&); \ INSTANTIATE_COMMS(TYPE) #ifdef CODI_FORWARD_TYPE diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 23bad36a8a4..9d19aa7c5a0 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1588,6 +1588,20 @@ void CFEASolver::BC_Clamped_Post(CGeometry *geometry, const CConfig *config, uns } +namespace { +/*--- Helper for BC_Sym_Plane ---*/ +template +void SubtractProjection(unsigned short nDim, const su2double* n, const Read& read, const Write& write) { + su2double tmp[3] = {}; + for (auto iDim = 0u; iDim < nDim; ++iDim) tmp[iDim] = read(iDim); + const su2double proj = DotProduct(nDim, tmp, n); + for (auto iDim = 0u; iDim < nDim; ++iDim) { + tmp[iDim] -= proj * n[iDim]; + write(iDim, tmp[iDim]); + } +} +} + void CFEASolver::BC_Sym_Plane(CGeometry *geometry, const CConfig *config, unsigned short val_marker) { if (geometry->GetnElem_Bound(val_marker) == 0) return; @@ -1610,14 +1624,8 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, const CConfig *config, unsign case 3: TriangleNormal(nodeCoord, normal); break; case 4: QuadrilateralNormal(nodeCoord, normal); break; } - - auto axis = 0u; - for (auto iDim = 1u; iDim < MAXNDIM; ++iDim) - axis = (fabs(normal[iDim]) > fabs(normal[axis]))? iDim : axis; - - if (fabs(normal[axis]) < 0.99*Norm(int(MAXNDIM),normal)) { - SU2_MPI::Error("The structural solver only supports axis-aligned symmetry planes.",CURRENT_FUNCTION); - } + const su2double area = Norm(int(MAXNDIM), normal); + for (auto iDim = 0u; iDim < MAXNDIM; ++iDim) normal[iDim] /= area; /*--- Impose zero displacement perpendicular to the symmetry plane. ---*/ @@ -1627,20 +1635,29 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, const CConfig *config, unsign const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Set and enforce solution at current and previous time-step ---*/ - nodes->SetSolution(iPoint, axis, 0.0); +#define SUBTRACT_PROJECTION(READ, WRITE) \ + SubtractProjection(nDim, normal, [&](unsigned short iDim) { return nodes->READ(iPoint, iDim); }, \ + [&](unsigned short iDim, const su2double& x) { nodes->WRITE(iPoint, iDim, x); }); + SUBTRACT_PROJECTION(GetSolution, SetSolution) if (dynamic) { - nodes->SetSolution_Vel(iPoint, axis, 0.0); - nodes->SetSolution_Accel(iPoint, axis, 0.0); - nodes->Set_Solution_time_n(iPoint, axis, 0.0); - nodes->SetSolution_Vel_time_n(iPoint, axis, 0.0); - nodes->SetSolution_Accel_time_n(iPoint, axis, 0.0); + SUBTRACT_PROJECTION(GetSolution_Vel, SetSolution_Vel) + SUBTRACT_PROJECTION(GetSolution_Accel, SetSolution_Accel) + SUBTRACT_PROJECTION(GetSolution_time_n, Set_Solution_time_n) + SUBTRACT_PROJECTION(GetSolution_Vel_time_n, SetSolution_Vel_time_n) + SUBTRACT_PROJECTION(GetSolution_Accel_time_n, SetSolution_Accel_time_n) } /*--- Set and enforce 0 solution for mesh deformation ---*/ - nodes->SetBound_Disp(iPoint, axis, 0.0); - LinSysSol(iPoint, axis) = 0.0; - if (LinSysReact.GetLocSize() > 0) LinSysReact(iPoint, axis) = 0.0; - Jacobian.EnforceSolutionAtDOF(iPoint, axis, su2double(0.0), LinSysRes); + SUBTRACT_PROJECTION(GetBound_Disp, SetBound_Disp) +#undef SUBTRACT_PROJECTION + + SubtractProjection(nDim, normal, [&](unsigned short iDim) { return LinSysSol(iPoint, iDim); }, + [&](unsigned short iDim, const su2double& x) { LinSysSol(iPoint, iDim) = x; }); + if (LinSysReact.GetLocSize() > 0) { + SubtractProjection(nDim, normal, [&](unsigned short iDim) { return LinSysReact(iPoint, iDim); }, + [&](unsigned short iDim, const su2double& x) { LinSysReact(iPoint, iDim) = x; }); + } + Jacobian.EnforceZeroProjection(iPoint, normal, LinSysRes); } diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 8bc45cb6a8b..530328b52a6 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1233,7 +1233,7 @@ def main(): # For a thin disk with the inner and outer radius of this geometry, from # "Formulas for Stress, Strain, and Structural Matrices", 2nd Edition, figure 19-4, # the maximum stress is 165.6MPa, we get a Von Misses stress very close to that. - rotating_cylinder_fea.test_vals = [-6.005467, -5.615543, -5.615527, 38, -8.126591, 1.6457e8] + rotating_cylinder_fea.test_vals = [-6.861940, -6.835545, -6.895500, 22, -8.313847, 1.6502e+08] test_list.append(rotating_cylinder_fea) # Dynamic beam, 2d