From 5fa37fb41cda43c8ea78c267b5bf5b41068aafc6 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 19 Dec 2022 11:47:15 +0000
Subject: [PATCH] BUG: fvMatrix: avoiding dupolicate adjustment. Fixes #2658

---
 src/finiteArea/faMatrices/faMatrix/faMatrix.C | 24 ++++++++++++-------
 .../fvMatrices/fvMatrix/fvMatrix.C            | 23 ++++++++++++------
 2 files changed, 32 insertions(+), 15 deletions(-)

diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.C b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
index b98f60e66a2..b2d81697989 100644
--- a/src/finiteArea/faMatrices/faMatrix/faMatrix.C
+++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
@@ -328,14 +328,11 @@ void Foam::faMatrix<Type>::setValuesFromList
             GeometricField<Type, faPatchField, areaMesh>&
         >(psi_).primitiveFieldRef();
 
-    forAll(faceLabels, i)
-    {
-        const label facei = faceLabels[i];
-        const Type& value = values[i];
-
-        psi[facei] = value;
-        source_[facei] = value*Diag[facei];
-    }
+    // Following actions:
+    // - adjust local field psi
+    // - set local matrix to be diagonal (so adjust source)
+    //      - cut connections to neighbours
+    // - make (on non-adjusted cells) contribution explicit
 
     if (symmetric() || asymmetric())
     {
@@ -392,6 +389,17 @@ void Foam::faMatrix<Type>::setValuesFromList
             }
         }
     }
+
+    // Note: above loop might have affected source terms on adjusted cells
+    // so make sure to adjust them afterwards
+    forAll(faceLabels, i)
+    {
+        const label facei = faceLabels[i];
+        const Type& value = values[i];
+
+        psi[facei] = value;
+        source_[facei] = value*Diag[facei];
+    }
 }
 
 
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 0a0879db777..a82341e6ce5 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -242,14 +242,12 @@ void Foam::fvMatrix<Type>::setValuesFromList
             GeometricField<Type, fvPatchField, volMesh>&
         >(psi_).primitiveFieldRef();
 
-    forAll(cellLabels, i)
-    {
-        const label celli = cellLabels[i];
-        const Type& value = values[i];
 
-        psi[celli] = value;
-        source_[celli] = value*Diag[celli];
-    }
+    // Following actions:
+    // - adjust local field psi
+    // - set local matrix to be diagonal (so adjust source)
+    //      - cut connections to neighbours
+    // - make (on non-adjusted cells) contribution explicit
 
     if (symmetric() || asymmetric())
     {
@@ -306,6 +304,17 @@ void Foam::fvMatrix<Type>::setValuesFromList
             }
         }
     }
+
+    // Note: above loop might have affected source terms on adjusted cells
+    // so make sure to adjust them afterwards
+    forAll(cellLabels, i)
+    {
+        const label celli = cellLabels[i];
+        const Type& value = values[i];
+
+        psi[celli] = value;
+        source_[celli] = value*Diag[celli];
+    }
 }
 
 
-- 
GitLab