Update to coupled solvers

applications/solvers/coupled/pUCoupledFoam/couplingTerms.H

@@ -0,0 +1,15 @@
+{
+    // Calculate grad p coupling matrix. Needs to be here if one uses
+    // gradient schemes with limiters.  VV, 9/June/2014
+    BlockLduSystem<vector, vector> pInU(fvm::grad(p));
+
+    // Calculate div U coupling.  Could be calculated only once since
+    // it is only geometry dependent.  VV, 9/June/2014
+    BlockLduSystem<vector, scalar> UInp(fvm::UDiv(U));
+
+    // Last argument in insertBlockCoupling says if the column direction
+    // should be incremented. This is needed for arbitrary positioning
+    // of U and p in the system. This could be better. VV, 30/April/2014
+    UpEqn.insertBlockCoupling(0, 3, pInU, true);
+    UpEqn.insertBlockCoupling(3, 0, UInp, false);
+}

applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C

@@ -38,9 +38,9 @@ Authors
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "fvBlockMatrix.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
-#include "fvBlockMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,9 +54,6 @@ int main(int argc, char *argv[])
 #   include "initContinuityErrs.H"
 #   include "readBlockSolverControls.H"
 
-    // Calculate div U coupling only once since it is only geometry dependant.
-    BlockLduSystem<vector, scalar> UInp(fvm::div(U));
-
     Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
@@ -73,15 +70,8 @@ int main(int argc, char *argv[])
         // Assemble and insert pressure equation
 #       include "pEqn.H"
 
-        // Calculate grad p coupling matrix. Needs to be here if one uses
+        // Assemble and insert coupling terms
-        // gradient schemes with limiters. VV, 9/June/2014
+#       include "couplingTerms.H"
-        BlockLduSystem<vector, vector> pInU(fvm::grad(p));
-
-        // Last argument in insertBlockCoupling says if the column direction
-        // should be incremented. This is needed for arbitrary positioning
-        // of U and p in the system. This could be better. VV, 30/April/2014
-        UpEqn.insertBlockCoupling(0, 3, pInU, true);
-        UpEqn.insertBlockCoupling(3, 0, UInp, false);
 
         // Solve the block matrix
         UpEqn.solve();

src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C

@@ -103,7 +103,7 @@ tmp
 <
     BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
 >
-divScheme<Type>::fvmDiv
+divScheme<Type>::fvmUDiv
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
@@ -128,6 +128,38 @@ divScheme<Type>::fvmDiv
 }
 
 
+template<class Type>
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+>
+divScheme<Type>::fvmUDiv
+(
+    const surfaceScalarField& flux,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    FatalErrorIn
+    (
+        "tmp<BlockLduSystem> divScheme<Type>::fvmDiv\n"
+        "(\n"
+        "    surfaceScalarField&"
+        "    GeometricField<Type, fvPatchField, volMesh>&"
+        ")\n"
+    )   << "Implicit div operator currently defined only for Gauss linear. "
+        << abort(FatalError);
+
+    typedef typename innerProduct<vector, Type>::type DivType;
+
+    tmp<BlockLduSystem<vector, DivType> > tbs
+    (
+        new BlockLduSystem<vector, DivType>(vf.mesh())
+    );
+
+    return tbs;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv

src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H

@@ -157,10 +157,19 @@ public:
         virtual tmp
         <
             BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
-        > fvmDiv
+        > fvmUDiv
         (
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const;
+
+        virtual tmp
+        <
+            BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+        > fvmUDiv
+        (
+            const surfaceScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
 };
 
 

src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C

@@ -72,14 +72,14 @@ template<class Type>
 tmp
 <
     BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
-> gaussDivScheme<Type>::fvmDiv
+> gaussDivScheme<Type>::fvmUDiv
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
    FatalErrorIn
    (
-       "tmp<BlockLduSystem> fvmDiv\n"
+       "tmp<BlockLduSystem> fvmUDiv\n"
        "(\n"
        "    GeometricField<Type, fvPatchField, volMesh>&"
        ")\n"
@@ -97,6 +97,37 @@ tmp
 }
 
 
+template<class Type>
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+> gaussDivScheme<Type>::fvmUDiv
+(
+    const surfaceScalarField& flux,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+   FatalErrorIn
+   (
+       "tmp<BlockLduSystem> fvmUDiv\n"
+       "(\n"
+       "    const surfaceScalarField& flux"
+       "    const GeometricField<Type, fvPatchField, volMesh>&"
+       ")\n"
+   )   << "Implicit div operator defined only for vector."
+       << abort(FatalError);
+
+   typedef typename innerProduct<vector, Type>::type DivType;
+
+   tmp<BlockLduSystem<vector, DivType> > tbs
+   (
+       new BlockLduSystem<vector, DivType>(vf.mesh())
+   );
+
+   return tbs;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv

src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H

@@ -102,10 +102,19 @@ public:
         tmp
         <
             BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
-        > fvmDiv
+        > fvmUDiv
         (
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const;
+
+        tmp
+        <
+            BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+        > fvmUDiv
+        (
+            const surfaceScalarField& flux,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
 };
 
 

src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C

@@ -42,7 +42,7 @@ namespace fv
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<>
-tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
+tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
 (
     const GeometricField<vector, fvPatchField, volMesh>& vf
 ) const
@@ -60,9 +60,9 @@ tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
     scalarField& source = bs.source();
 
     // Grab ldu parts of block matrix as linear always
-    typename CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
+    CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
-    typename CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
+    CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
-    typename CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
+    CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
 
     const vectorField& SfIn = mesh.Sf().internalField();
 
@@ -89,9 +89,9 @@ tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
 
         if (patch.coupled())
         {
-            typename CoeffField<vector>::linearTypeField& pcoupleUpper =
+            CoeffField<vector>::linearTypeField& pcoupleUpper =
                 bs.coupleUpper()[patchI].asLinear();
-            typename CoeffField<vector>::linearTypeField& pcoupleLower =
+            CoeffField<vector>::linearTypeField& pcoupleLower =
                 bs.coupleLower()[patchI].asLinear();
 
             const vectorField pcl = -pw*Sf;
@@ -119,6 +119,95 @@ tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
 }
 
 
+template<>
+tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
+(
+    const surfaceScalarField& flux,
+    const GeometricField<vector, fvPatchField, volMesh>& vf
+) const
+{
+    tmp<surfaceScalarField> tweights = this->tinterpScheme_().weights(vf);
+    const scalarField& wIn = tweights().internalField();
+
+    const fvMesh& mesh = vf.mesh();
+
+    tmp<BlockLduSystem<vector, scalar> > tbs
+    (
+        new BlockLduSystem<vector, scalar>(mesh)
+    );
+    BlockLduSystem<vector, scalar>& bs = tbs();
+    scalarField& source = bs.source();
+
+    // Grab ldu parts of block matrix as linear always
+    CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
+    CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
+    CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
+
+    const vectorField& SfIn = mesh.Sf().internalField();
+
+    const scalarField& fluxIn = flux.internalField();
+
+    l = -wIn*fluxIn*SfIn;
+    u = l + fluxIn*SfIn;
+    bs.negSumDiag();
+
+    // Boundary contributions
+    forAll(vf.boundaryField(), patchI)
+    {
+        const fvPatchVectorField& pf = vf.boundaryField()[patchI];
+        const fvPatch& patch = pf.patch();
+        const vectorField& Sf = patch.Sf();
+        const fvsPatchScalarField& pw = tweights().boundaryField()[patchI];
+        const unallocLabelList& fc = patch.faceCells();
+
+        const scalarField& pFlux = flux.boundaryField()[patchI];
+
+        const vectorField internalCoeffs(pf.valueInternalCoeffs(pw));
+
+        // Diag contribution
+        forAll(pf, faceI)
+        {
+            d[fc[faceI]] += cmptMultiply
+            (
+                internalCoeffs[faceI],
+                pFlux[faceI]*Sf[faceI]
+            );
+        }
+
+        if (patch.coupled())
+        {
+            CoeffField<vector>::linearTypeField& pcoupleUpper =
+                bs.coupleUpper()[patchI].asLinear();
+            CoeffField<vector>::linearTypeField& pcoupleLower =
+                bs.coupleLower()[patchI].asLinear();
+
+            const vectorField pcl = -pw*pFlux*Sf;
+            const vectorField pcu = pcl + pFlux*Sf;
+
+            // Coupling  contributions
+            pcoupleLower -= pcl;
+            pcoupleUpper -= pcu;
+        }
+        else
+        {
+            const vectorField boundaryCoeffs(pf.valueBoundaryCoeffs(pw));
+
+            // Boundary contribution
+            forAll(pf, faceI)
+            {
+                source[fc[faceI]] -=
+                (
+                    boundaryCoeffs[faceI] & (pFlux[faceI]*Sf[faceI])
+                );
+            }
+        }
+    }
+
+    // Interpolation schemes with corrections not accounted for
+
+    return tbs;
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv

src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H

@@ -51,11 +51,18 @@ namespace fv
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<>
-tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
+tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
 (
     const GeometricField<vector, fvPatchField, volMesh>&
 ) const;
 
+template<>
+tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
+(
+    const surfaceScalarField& flux,
+    const GeometricField<vector, fvPatchField, volMesh>&
+) const;
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteVolume/finiteVolume/fvm/fvmDiv.C

@@ -103,7 +103,7 @@ template<class Type>
 tmp
 <
     BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
-> div
+> UDiv
 (
     GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
@@ -113,7 +113,7 @@ tmp
     (
         vf.mesh(),
         vf.mesh().schemesDict().divScheme(name)
-    )().fvmDiv(vf);
+    )().fvmUDiv(vf);
 }
 
 
@@ -121,12 +121,52 @@ template<class Type>
 tmp
 <
     BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
-> div
+> UDiv
+(
+    const surfaceScalarField& flux,
+    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const word& name
+)
+{
+    return fv::divScheme<Type>::New
+    (
+        vf.mesh(),
+        vf.mesh().schemesDict().divScheme(name)
+    )().fvmUDiv(flux, vf);
+}
+
+
+template<class Type>
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+> UDiv
+(
+    const tmp<surfaceScalarField>& tflux,
+    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const word& name
+)
+{
+    tmp
+    <
+        BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+    >
+    Div(fvm::UDiv(tflux(), vf, name));
+    tflux.clear();
+    return Div;
+}
+
+
+template<class Type>
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+> UDiv
 (
     GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
-    return fvm::div
+    return fvm::UDiv
     (
         vf,
         "div(" + vf.name() + ')'
@@ -134,6 +174,45 @@ tmp
 }
 
 
+template<class Type>
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+> UDiv
+(
+    const surfaceScalarField& flux,
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    return fvm::UDiv
+    (
+        flux,
+        vf,
+        "div(" + vf.name() + ')'
+    );
+}
+
+
+template<class Type>
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+> UDiv
+(
+    const tmp<surfaceScalarField>& tflux,
+    GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp
+    <
+        BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+    >
+    Div(fvm::UDiv(tflux(), vf));
+    tflux.clear();
+    return Div;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fvm

src/finiteVolume/finiteVolume/fvm/fvmDiv.H

@@ -68,7 +68,6 @@ namespace fvm
         const word& name
     );
 
-
     template<class Type>
     tmp<fvMatrix<Type> > div
     (
@@ -83,12 +82,12 @@ namespace fvm
         GeometricField<Type, fvPatchField, volMesh>&
     );
 
-    // Implicit div operator for block systems
+    // Implicit div operators for block systems
     template<class Type>
     tmp
     <
         BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
-    > div
+    > UDiv
     (
         GeometricField<Type, fvPatchField, volMesh>&,
         const word&
@@ -98,8 +97,50 @@ namespace fvm
     tmp
     <
         BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
-    > div
+    > UDiv
     (
+        const surfaceScalarField&,
+        GeometricField<Type, fvPatchField, volMesh>&,
+        const word&
+    );
+
+    template<class Type>
+    tmp
+    <
+        BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+    > UDiv
+    (
+        const tmp<surfaceScalarField>&,
+        GeometricField<Type, fvPatchField, volMesh>&,
+        const word&
+    );
+
+    template<class Type>
+    tmp
+    <
+        BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+    > UDiv
+    (
+        GeometricField<Type, fvPatchField, volMesh>&
+    );
+
+    template<class Type>
+    tmp
+    <
+        BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+    > UDiv
+    (
+        const surfaceScalarField&,
+        GeometricField<Type, fvPatchField, volMesh>&
+    );
+
+    template<class Type>
+    tmp
+    <
+        BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+    > UDiv
+    (
+        const tmp<surfaceScalarField>&,
         GeometricField<Type, fvPatchField, volMesh>&
     );
 }

src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.C

@@ -60,9 +60,9 @@ tmp<BlockLduSystem<vector, vector> > gaussGrad<scalar>::fvmGrad
     vectorField& source = bs.source();
 
     // Grab ldu parts of block matrix as linear always
-    typename CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
+    CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
-    typename CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
+    CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
-    typename CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
+    CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
 
     const vectorField& SfIn = mesh.Sf().internalField();
 
@@ -89,9 +89,9 @@ tmp<BlockLduSystem<vector, vector> > gaussGrad<scalar>::fvmGrad
 
         if (patch.coupled())
         {
-            typename CoeffField<vector>::linearTypeField& pcoupleUpper =
+            CoeffField<vector>::linearTypeField& pcoupleUpper =
                 bs.coupleUpper()[patchI].asLinear();
-            typename CoeffField<vector>::linearTypeField& pcoupleLower =
+            CoeffField<vector>::linearTypeField& pcoupleLower =
                 bs.coupleLower()[patchI].asLinear();
 
             const vectorField pcl = -pw*Sf;

src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C

@@ -80,9 +80,9 @@ tmp<BlockLduSystem<vector, vector> > leastSquaresGrad<scalar>::fvmGrad
     vectorField& source = bs.source();
 
     // Grab ldu parts of block matrix as linear always
-    typename CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
+    CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
-    typename CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
+    CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
-    typename CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
+    CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
 
     // Get references to least square vectors
     const leastSquaresVectors& lsv = leastSquaresVectors::New(mesh);
@@ -125,9 +125,9 @@ tmp<BlockLduSystem<vector, vector> > leastSquaresGrad<scalar>::fvmGrad
             const scalarField cellVInNei =
                 cellV.boundaryField()[patchI].patchNeighbourField();
 
-            typename CoeffField<vector>::linearTypeField& pcoupleUpper =
+            CoeffField<vector>::linearTypeField& pcoupleUpper =
                 bs.coupleUpper()[patchI].asLinear();
-            typename CoeffField<vector>::linearTypeField& pcoupleLower =
+            CoeffField<vector>::linearTypeField& pcoupleLower =
                 bs.coupleLower()[patchI].asLinear();
 
             // Coupling  and diagonal contributions