Temporary implicit block boundary conditions and porous/MRF

applications/solvers/coupled/MRFPorousFoam/UEqn.H

@@ -53,80 +53,6 @@
     // Insert momentum equation
     UpEqn.insertEquation(0, UEqn);
 
-    if (addMRF || addPorosity)
+#   include "addImplicitMRFPorous.H"
-    {
+#   include "addBlockCoupledBC.H"
-        // Manually over-ride the 3x3 block to handle the off-diagonal
-        // part of the Ap coefficient
-        const tensorField& TUIn = tTU().internalField();
-
-        CoeffField<vector4>::squareTypeField& DD = UpEqn.diag().asSquare();
-
-        const scalarField& V = mesh.V().field();
-
-        // Note: insertion should only happen in MRF and porous cell zones
-        // HJ, 14/Mar/2016
-
-        // Warning.  Possible problem with a zone which is both MRF and porous
-        // The solution is to do the loop below everywhere
-        // Currently only inserting zones for efficiency.  HJ, 18/Nov/2017
-        register label cellI;
-
-        forAll (mrfZones, mrfZoneI)
-        {
-            const labelList& curZoneCells = mrfZones[mrfZoneI].zone();
-
-            // Loop over all cells in the zone
-            forAll (curZoneCells, zcI)
-            {
-                cellI = curZoneCells[zcI];
-
-                const scalar& cellV = V[cellI];
-
-                const tensor& cellTU = TUIn[cellI];
-
-                CoeffField<vector4>::squareType& cellDD = DD[cellI];
-
-                cellDD(0, 0) += cellV*cellTU.xx();
-                cellDD(0, 1) += cellV*cellTU.xy();
-                cellDD(0, 2) += cellV*cellTU.xz();
-
-                cellDD(1, 0) += cellV*cellTU.yx();
-                cellDD(1, 1) += cellV*cellTU.yy();
-                cellDD(2, 2) += cellV*cellTU.yz();
-
-                cellDD(2, 0) += cellV*cellTU.zx();
-                cellDD(2, 1) += cellV*cellTU.zy();
-                cellDD(2, 2) += cellV*cellTU.zz();
-            }
-        }
-
-        forAll (pZones, pZoneI)
-        {
-            const labelList& curZoneCells = pZones[pZoneI].zone();
-
-            // Loop over all cells in the zone
-            forAll (curZoneCells, zcI)
-            {
-                cellI = curZoneCells[zcI];
-
-                const scalar& cellV = V[cellI];
-
-                const tensor& cellTU = TUIn[cellI];
-
-                CoeffField<vector4>::squareType& cellDD = DD[cellI];
-
-                cellDD(0, 0) += cellV*cellTU.xx();
-                cellDD(0, 1) += cellV*cellTU.xy();
-                cellDD(0, 2) += cellV*cellTU.xz();
-
-                cellDD(1, 0) += cellV*cellTU.yx();
-                cellDD(1, 1) += cellV*cellTU.yy();
-                cellDD(2, 2) += cellV*cellTU.yz();
-
-                cellDD(2, 0) += cellV*cellTU.zx();
-                cellDD(2, 1) += cellV*cellTU.zy();
-                cellDD(2, 2) += cellV*cellTU.zz();
-            }
-        }
-    }
 }

applications/solvers/coupled/MRFPorousFoam/addBlockCoupledBC.H

@@ -0,0 +1,56 @@
+    // Hack block-coupled boundary conditions: due for rewrite
+    const volScalarField nuEff = turbulence->nuEff();
+
+    forAll (U.boundaryField(), patchI)
+    {
+        if (U.boundaryField()[patchI].blockCoupled())
+        {
+            // Insert correcting fully implicit coupling coefficient
+
+            const labelList fc = mesh.boundary()[patchI].faceCells();
+
+            const fvPatchVectorField& Up = U.boundaryField()[patchI];
+
+            // Warning: hacked for nuEff in viscosity
+            const scalarField nutpMagSf =
+                nuEff.boundaryField()[patchI]*
+                mesh.magSf().boundaryField()[patchI];
+
+            // Get boundary condition contribution to matrix diagonal
+            tensorField patchDiag =
+                -Up.blockGradientInternalCoeffs()().asSquare()*nutpMagSf;
+
+            // Get matrix diagonal
+            CoeffField<vector4>::squareTypeField& blockDiag =
+                UpEqn.diag().asSquare();
+
+            forAll (fc, faceI)
+            {
+                blockDiag[fc[faceI]](0, 0) += patchDiag[faceI].xx();
+                blockDiag[fc[faceI]](0, 1) += patchDiag[faceI].xy();
+                blockDiag[fc[faceI]](0, 2) += patchDiag[faceI].xz();
+
+                blockDiag[fc[faceI]](1, 0) += patchDiag[faceI].yx();
+                blockDiag[fc[faceI]](1, 1) += patchDiag[faceI].yy();
+                blockDiag[fc[faceI]](1, 2) += patchDiag[faceI].yz();
+
+                blockDiag[fc[faceI]](2, 0) += patchDiag[faceI].zx();
+                blockDiag[fc[faceI]](3, 1) += patchDiag[faceI].zy();
+                blockDiag[fc[faceI]](3, 2) += patchDiag[faceI].zz();
+            }
+
+         // Get boundary condition contribution to matrix source
+         vectorField patchSource =
+             -Up.blockGradientBoundaryCoeffs()*nutpMagSf;
+
+         // Get matrix source
+         Field<vector4>& blockSource = UpEqn.source();
+
+            forAll (fc, faceI)
+            {
+                blockSource[fc[faceI]](0) -= patchSource[faceI](0);
+                blockSource[fc[faceI]](1) -= patchSource[faceI](1);
+                blockSource[fc[faceI]](2) -= patchSource[faceI](2);
+            }
+        }
+    }

applications/solvers/coupled/MRFPorousFoam/addImplicitMRFPorous.H

@@ -0,0 +1,77 @@
+
+    if (addMRF || addPorosity)
+    {
+        // Manually over-ride the 3x3 block to handle the off-diagonal
+        // part of the Ap coefficient
+        const tensorField& TUIn = tTU().internalField();
+
+        CoeffField<vector4>::squareTypeField& DD = UpEqn.diag().asSquare();
+
+        const scalarField& V = mesh.V().field();
+
+        // Note: insertion should only happen in MRF and porous cell zones
+        // HJ, 14/Mar/2016
+
+        // Warning.  Possible problem with a zone which is both MRF and porous
+        // The solution is to do the loop below everywhere
+        // Currently only inserting zones for efficiency.  HJ, 18/Nov/2017
+        register label cellI;
+
+        forAll (mrfZones, mrfZoneI)
+        {
+            const labelList& curZoneCells = mrfZones[mrfZoneI].zone();
+
+            // Loop over all cells in the zone
+            forAll (curZoneCells, zcI)
+            {
+                cellI = curZoneCells[zcI];
+
+                const scalar& cellV = V[cellI];
+
+                const tensor& cellTU = TUIn[cellI];
+
+                CoeffField<vector4>::squareType& cellDD = DD[cellI];
+
+                cellDD(0, 0) += cellV*cellTU.xx();
+                cellDD(0, 1) += cellV*cellTU.xy();
+                cellDD(0, 2) += cellV*cellTU.xz();
+
+                cellDD(1, 0) += cellV*cellTU.yx();
+                cellDD(1, 1) += cellV*cellTU.yy();
+                cellDD(2, 2) += cellV*cellTU.yz();
+
+                cellDD(2, 0) += cellV*cellTU.zx();
+                cellDD(2, 1) += cellV*cellTU.zy();
+                cellDD(2, 2) += cellV*cellTU.zz();
+            }
+        }
+
+        forAll (pZones, pZoneI)
+        {
+            const labelList& curZoneCells = pZones[pZoneI].zone();
+
+            // Loop over all cells in the zone
+            forAll (curZoneCells, zcI)
+            {
+                cellI = curZoneCells[zcI];
+
+                const scalar& cellV = V[cellI];
+
+                const tensor& cellTU = TUIn[cellI];
+
+                CoeffField<vector4>::squareType& cellDD = DD[cellI];
+
+                cellDD(0, 0) += cellV*cellTU.xx();
+                cellDD(0, 1) += cellV*cellTU.xy();
+                cellDD(0, 2) += cellV*cellTU.xz();
+
+                cellDD(1, 0) += cellV*cellTU.yx();
+                cellDD(1, 1) += cellV*cellTU.yy();
+                cellDD(2, 2) += cellV*cellTU.yz();
+
+                cellDD(2, 0) += cellV*cellTU.zx();
+                cellDD(2, 1) += cellV*cellTU.zy();
+                cellDD(2, 2) += cellV*cellTU.zz();
+            }
+        }
+    }

applications/solvers/coupled/pUCoupledFoam/UEqn.H

@@ -12,4 +12,6 @@
     UEqn.relax();
 
     UpEqn.insertEquation(0, UEqn);
+
+#   include "addBlockCoupledBC.H"
 }

applications/solvers/coupled/pUCoupledFoam/addBlockCoupledBC.H

@@ -0,0 +1,56 @@
+    // Hack block-coupled boundary conditions: due for rewrite
+    const volScalarField nuEff = turbulence->nuEff();
+
+    forAll (U.boundaryField(), patchI)
+    {
+        if (U.boundaryField()[patchI].blockCoupled())
+        {
+            // Insert correcting fully implicit coupling coefficient
+
+            const labelList fc = mesh.boundary()[patchI].faceCells();
+
+            const fvPatchVectorField& Up = U.boundaryField()[patchI];
+
+            // Warning: hacked for nuEff in viscosity
+            const scalarField nutpMagSf =
+                nuEff.boundaryField()[patchI]*
+                mesh.magSf().boundaryField()[patchI];
+
+            // Get boundary condition contribution to matrix diagonal
+            tensorField patchDiag =
+                -Up.blockGradientInternalCoeffs()().asSquare()*nutpMagSf;
+
+            // Get matrix diagonal
+            CoeffField<vector4>::squareTypeField& blockDiag =
+                UpEqn.diag().asSquare();
+
+            forAll (fc, faceI)
+            {
+                blockDiag[fc[faceI]](0, 0) += patchDiag[faceI].xx();
+                blockDiag[fc[faceI]](0, 1) += patchDiag[faceI].xy();
+                blockDiag[fc[faceI]](0, 2) += patchDiag[faceI].xz();
+
+                blockDiag[fc[faceI]](1, 0) += patchDiag[faceI].yx();
+                blockDiag[fc[faceI]](1, 1) += patchDiag[faceI].yy();
+                blockDiag[fc[faceI]](1, 2) += patchDiag[faceI].yz();
+
+                blockDiag[fc[faceI]](2, 0) += patchDiag[faceI].zx();
+                blockDiag[fc[faceI]](3, 1) += patchDiag[faceI].zy();
+                blockDiag[fc[faceI]](3, 2) += patchDiag[faceI].zz();
+            }
+
+         // Get boundary condition contribution to matrix source
+         vectorField patchSource =
+             -Up.blockGradientBoundaryCoeffs()*nutpMagSf;
+
+         // Get matrix source
+         Field<vector4>& blockSource = UpEqn.source();
+
+            forAll (fc, faceI)
+            {
+                blockSource[fc[faceI]](0) -= patchSource[faceI](0);
+                blockSource[fc[faceI]](1) -= patchSource[faceI](1);
+                blockSource[fc[faceI]](2) -= patchSource[faceI](2);
+            }
+        }
+    }