rhoSimplecFoam: Added coupled-patch handling to H1 to support parallel running

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -54,7 +54,8 @@ else
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rho*rAU, p) == fvc::div(phiHbyA)
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rho*rAU, p)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -63,7 +64,7 @@ else
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA - pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -88,5 +89,10 @@ if (closedVolume)
 rho = thermo.rho();
 rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
-rho.relax();
+
+if (!simple.transonic())
+{
+    rho.relax();
+}
+
 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H

@@ -3,8 +3,6 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField p0(p);
-
 volScalarField AU(UEqn().A());
 volScalarField AtU(AU - UEqn().H1());
 
@@ -22,7 +20,7 @@ if (simple.transonic())
         surfaceScalarField phid
         (
             "phid",
-            fvc::interpolate(psi*HbyA) & mesh.Sf()
+            fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
         );
 
         surfaceScalarField phic
@@ -31,6 +29,8 @@ if (simple.transonic())
             fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf()
         );
 
+        HbyA -= fvc::grad(p)*(1.0/AU - 1.0/AtU);
+
         fvScalarMatrix pEqn
         (
             fvm::div(phid, p)
@@ -58,11 +58,15 @@ else
         surfaceScalarField phiHbyA
         (
             "phiHbyA",
-            fvc::interpolate(rho*HbyA) & mesh.Sf()
+            fvc::interpolate(rho)*fvc::interpolate(HbyA) & mesh.Sf()
         );
 
-        closedVolume = adjustPhi(phi, U, p);
-        phi += fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf();
+        closedVolume = adjustPhi(phiHbyA, U, p);
+
+        phiHbyA +=
+            fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf();
+
+        HbyA -= fvc::grad(p)*(1.0/AU - 1.0/AtU);
 
         fvScalarMatrix pEqn
         (
@@ -88,9 +92,7 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U = HbyA - (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU);
-//U = HbyA - fvc::grad(p)/AU;
-
+U = HbyA - fvc::grad(p)/AtU;
 U.correctBoundaryConditions();
 
 // For closed-volume cases adjust the pressure and density levels

src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C

@@ -99,8 +99,16 @@ void Foam::SolverPerformance<Type>::print
 {
     for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
     {
-        os  << solverName_ << ":  Solving for "
-            << word(fieldName_ + pTraits<Type>::componentNames[cmpt]);
+        if (pTraits<Type>::nComponents == 1)
+        {
+            os  << solverName_ << ":  Solving for " << fieldName_;
+
+        }
+        else
+        {
+            os  << solverName_ << ":  Solving for "
+                << word(fieldName_ + pTraits<Type>::componentNames[cmpt]);
+        }
 
         if (singular_[cmpt])
         {

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C

@@ -295,4 +295,36 @@ Foam::tmp<Foam::scalarField> Foam::lduMatrix::residual
 }
 
 
+Foam::tmp<Foam::scalarField > Foam::lduMatrix::H1() const
+{
+    tmp<scalarField > tH1
+    (
+        new scalarField(lduAddr().size(), 0.0)
+    );
+
+    if (lowerPtr_ || upperPtr_)
+    {
+        scalarField& H1_ = tH1();
+
+        scalar* __restrict__ H1Ptr = H1_.begin();
+
+        const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
+        const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
+
+        const scalar* __restrict__ lowerPtr = lower().begin();
+        const scalar* __restrict__ upperPtr = upper().begin();
+
+        register const label nFaces = upper().size();
+
+        for (register label face=0; face<nFaces; face++)
+        {
+            H1Ptr[uPtr[face]] -= lowerPtr[face];
+            H1Ptr[lPtr[face]] -= upperPtr[face];
+        }
+    }
+
+    return tH1;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C

@@ -351,36 +351,4 @@ void Foam::lduMatrix::operator*=(scalar s)
 }
 
 
-Foam::tmp<Foam::scalarField > Foam::lduMatrix::H1() const
-{
-    tmp<scalarField > tH1
-    (
-        new scalarField(lduAddr().size(), 0.0)
-    );
-
-    if (lowerPtr_ || upperPtr_)
-    {
-        scalarField& H1_ = tH1();
-
-        scalar* __restrict__ H1Ptr = H1_.begin();
-
-        const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
-        const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
-
-        const scalar* __restrict__ lowerPtr = lower().begin();
-        const scalar* __restrict__ upperPtr = upper().begin();
-
-        register const label nFaces = upper().size();
-
-        for (register label face=0; face<nFaces; face++)
-        {
-            H1Ptr[uPtr[face]] -= lowerPtr[face];
-            H1Ptr[lPtr[face]] -= upperPtr[face];
-        }
-    }
-
-    return tH1;
-}
-
-
 // ************************************************************************* //

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -841,25 +841,23 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::H1() const
     );
     volScalarField& H1_ = tH1();
 
-    // Loop over field components
-    /*
-    for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
-    {
-        scalarField psiCmpt(psi_.internalField().component(cmpt));
+    H1_.internalField() = lduMatrix::H1();
 
-        scalarField boundaryDiagCmpt(psi_.size(), 0.0);
-        addBoundaryDiag(boundaryDiagCmpt, cmpt);
-        boundaryDiagCmpt.negate();
-        addCmptAvBoundaryDiag(boundaryDiagCmpt);
+    forAll(psi_.boundaryField(), patchI)
+    {
+        const fvPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 
-        H1_.internalField().replace(cmpt, boundaryDiagCmpt);
+        if (ptf.coupled() && ptf.size())
+        {
+            addToInternalField
+            (
+                lduAddr().patchAddr(patchI),
+                boundaryCoeffs_[patchI].component(0),
+                H1_
+            );
+        }
     }
 
-    H1_.internalField() += lduMatrix::H1();
-    */
-
-    H1_.internalField() = lduMatrix::H1();
-
     H1_.internalField() /= psi_.mesh().V();
     H1_.correctBoundaryConditions();
 

tutorials/compressible/rhoSimplecFoam/squareBend/0/U

@@ -28,7 +28,7 @@ boundaryField
     inlet
     {
         type            flowRateInletVelocity;
-        flowRate        constant 0.5; //0.75;
+        flowRate        constant 0.5;
         value           uniform (0 0 0);
     }
     outlet

tutorials/compressible/rhoSimplecFoam/squareBend/system/decomposeParDict

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          hierarchical;
+
+simpleCoeffs
+{
+    n               (8 1 1);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (4 2 1);
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //