Skip to content

Commit 9d235a4

Browse files
committed
manual asymmetric block update
1 parent 9ec0e83 commit 9d235a4

File tree

1 file changed

+28
-2
lines changed

1 file changed

+28
-2
lines changed

SU2_CFD/src/solvers/CTurbSASolver.cpp

Lines changed: 28 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -337,8 +337,21 @@ void CTurbSASolver::Viscous_Residual(const unsigned long iEdge, const CGeometry*
337337
// }
338338
// else {
339339
LinSysRes.SubtractBlock(iPoint, residual_ij);
340+
// LinSysRes.AddBlock(jPoint, residual_ij); **Is this needed?
340341
if (implicit) {
341-
Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual_ij.jacobian_i, residual_ij.jacobian_j);
342+
auto* Block_ii = Jacobian.GetBlock(iPoint, iPoint);
343+
344+
for (int iVar=0; iVar<nVar; iVar++)
345+
for (int jVar=0; jVar<nVar; jVar++) {
346+
Block_ii[iVar*nVar + jVar] -= SU2_TYPE::GetValue(residual_ij.jacobian_i[iVar][jVar]);
347+
}
348+
349+
auto* Block_ij = Jacobian.GetBlock(iPoint, jPoint);
350+
for (int iVar=0; iVar<nVar; iVar++)
351+
for (int jVar=0; jVar<nVar; jVar++) {
352+
Block_ij[iVar*nVar + jVar] -= SU2_TYPE::GetValue(residual_ij.jacobian_j[iVar][jVar]);
353+
}
354+
// Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual_ij.jacobian_i, residual_ij.jacobian_j);
342355
// }
343356
}
344357

@@ -356,8 +369,21 @@ void CTurbSASolver::Viscous_Residual(const unsigned long iEdge, const CGeometry*
356369
// }
357370
// else {
358371
LinSysRes.SubtractBlock(jPoint, residual_ji);
372+
// LinSysRes.AddBlock(iPoint, residual_ji); **Is this needed?
359373
if (implicit) {
360-
Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual_ji.jacobian_i, residual_ji.jacobian_j);
374+
auto* Block_ji = Jacobian.GetBlock(jPoint, iPoint);
375+
//the order of arguments were flipped in the evaluation of residual_ji, the jacobian associated with point i is stored in jacobian_j and point j in jacobian_i
376+
for (int iVar=0; iVar<nVar; iVar++)
377+
for (int jVar=0; jVar<nVar; jVar++) {
378+
Block_ji[iVar*nVar + jVar] -= SU2_TYPE::GetValue(residual_ij.jacobian_j[iVar][jVar]);
379+
}
380+
381+
auto* Block_jj = Jacobian.GetBlock(jPoint, jPoint);
382+
for (int iVar=0; iVar<nVar; iVar++)
383+
for (int jVar=0; jVar<nVar; jVar++) {
384+
Block_jj[iVar*nVar + jVar] -= SU2_TYPE::GetValue(residual_ij.jacobian_i[iVar][jVar]);
385+
}
386+
// Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual_ji.jacobian_i, residual_ji.jacobian_j);
361387
// }
362388
}
363389
}

0 commit comments

Comments
 (0)