Merge pull request #8 from NREL/master

Merge changes from master

Author: Bart Doekemeijer

applications/solvers/incompressible/windEnergy/ABLSolver/Make/options

@@ -6,8 +6,9 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
-    -I./interpolate2D \
-    -I./windRoseToCartesian
+    -I../commonAlgorithms \
+    -I../commonAlgorithms/interpolate2D \
+    -I../commonAlgorithms/windRoseToCartesian
 
 
 EXE_LIBS = \

applications/solvers/incompressible/windEnergy/ABLSolver/UEqn.H

@@ -2,6 +2,8 @@
 
     #include "computeCoriolisForce.H"
 
+    #include "computeBuoyancyTerm.H"
+
     fvVectorMatrix UEqn
     (
         fvm::ddt(U)                             // time derivative
@@ -24,7 +26,7 @@
             (
                 (
                   - fvc::snGrad(p_rgh)          // modified pressure gradient
-                  - ghf*fvc::snGrad(rhok)       // buoyancy force
+                  + buoyancyTerm                // buoyancy force
                 ) * mesh.magSf()
             )
         );

applications/solvers/incompressible/windEnergy/ABLSolver/computeDivergence.H

@@ -136,10 +136,6 @@
         mesh
     );
 
-    // Create and calculate the gravity potential field
-    Info << "Creating and calculating the gravity potential field, gh and ghf..." << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
 
     // Create and calculate the static pressure field
     Info << "Creating and calculating the static pressure field, p..." << endl;
@@ -153,9 +149,11 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        p_rgh + rhok*gh
+      //p_rgh + rhok*gh
+        p_rgh
     );
 
+
     // Set up the pressure reference cell information
     Info << "Setting up the pressure reference cell information..." << endl;
     label pRefCell = 0;
@@ -169,13 +167,21 @@
         pRefValue
     );
 
-    if (p_rgh.needReference())
+
+    // Create and calculate the gravity potential field
+    Info << "Creating and calculating the gravity potential field, gh and ghf..." << endl;
+    
+    vector hRef_ = vector::zero;
+    if (pRefCell != -1)
     {
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
-        );
+        hRef_ = mesh.C()[pRefCell];
     }
+    reduce(hRef_,sumOp<vector>());
+    dimensionedVector hRef("hRef",dimLength,hRef_);
+
+    volScalarField gh("gh", g & (mesh.C() - hRef));
+    surfaceScalarField ghf("ghf", g & (mesh.Cf() - hRef));
+
 
+    // Compute the full pressure and adjust the pressure level.
+    #include "adjustPressureLevel.H"

applications/solvers/incompressible/windEnergy/ABLSolver/createFields.H

@@ -1,22 +1,30 @@
 {
+    // The inverse of the momentum equation "A" operator (the matrix diagonal) and
+    // its projection from cell centers to cell faces.  Because it is a diagonal matrix,
+    // The inverse is just the diagonal matrix of the reciprocals, hence the "r" in the
+    // name "rAU".
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
 
+    // The momentum equation "H" operator is the off-diagonal times the last known U.
+    // In the equations it is always multiplied with inv(A).  This is HbyA.
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
+    // This is the Boussinesq buoyancy term multipled by the inverse of the A operator.
+    #include "computeBuoyancyTerm.H"
+    surfaceScalarField phig(rAUf * buoyancyTerm * mesh.magSf());
 
+    // Project HbyA to cell faces and apply a correction for time stepping scheme.
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
       + rAUf*fvc::ddtCorr(U, phi)
     );
 
-    adjustPhi(phiHbyA, U, p_rgh);
-
-    surfaceScalarField phiFixedFlux = fvc::interpolate(U) & mesh.Sf();
+    // This part balances global mass flux.
+    surfaceScalarField phiFixedFlux = phiHbyA;
 
     adjustPhi(phiFixedFlux, U, p_rgh);
 
@@ -32,14 +40,19 @@
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
     );
 
+
+    // Non-orthogonal corrector loop.
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
             fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        if (activatePressureRefCell)
+        {
+            p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        }
 
         p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
@@ -58,16 +71,6 @@
         }
     }
 
-    p = p_rgh + rhok*gh;
-
-    if (p_rgh.needReference())
-    {
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
-        );
-        p_rgh = p - rhok*gh;
-    }
+    // Compute the full pressure and adjust the pressure level.
+    #include "adjustPressureLevel.H"
 }

applications/solvers/incompressible/windEnergy/ABLSolver/pEqn.H

@@ -198,6 +198,55 @@
 
 
 
+
+    // PROPERTIES CONCERNING THE WAY IN WHICH PERTURBATION PRESSURE IS DEFINED
+
+       // Options for defining the perturbation pressure:
+       // - noSplit:   do not split out hydrostatic part; pressure is then perturbation pressure.
+       // - rho0Split: split out the hydrostatic part; define hydrostatic as rho_0 * g * z.
+       // - rhokSplit: split out the hydrostatic part; define hydrostatic as rho_k * g * z.
+       word perturbationPressureType(ABLProperties.lookupOrDefault<word>("perturbationPressureType","rhokSplit"));
+       word perturbationOutput;
+       if (perturbationPressureType == "noSplit")
+       {
+           perturbationOutput = "nothing";
+       }
+       else if (perturbationPressureType == "rho0Split")
+       {
+           perturbationOutput = "rho_0 * g * z";
+       }
+       else if (perturbationPressureType == "rhokSplit")
+       {
+           perturbationOutput = "rho_k * g * z";
+       }
+       Info << "Defining background hydrostatic pressure to be " << perturbationOutput << endl;
+
+
+       // This switch dictates whether or not the pressure reference cell is set explicitly
+       // in the p_rgh solve.  If true, pressure is constrained at the pressure reference
+       // cell by manipulating the matrix row for that cell to remove the null space.  This
+       // assures that the pressure level is constrained, but it also means the continuity
+       // equation is no longer solved at that cell, so the divergence error can be significant
+       // enough there to cause localized oscillations.  Iterative solvers can still converge
+       // even with a null space, but more iterations may be needed, so this switch can be 
+       // set to false.
+       bool activatePressureRefCell(ABLProperties.lookupOrDefault<bool>("activatePressureRefCell", true));
+       if (activatePressureRefCell)
+       {
+            Info << "Pressure reference cell matrix row modification enabled" << endl;
+       }
+       else
+       {
+            Info << "Pressure reference cell matrix row modification not enabled" << endl;
+       }
+
+
+
+     
+
+
+
+
     // PROPERTIES CONCERNING GATHERING STATISTICS
 
        // Gather/write statistics?

applications/solvers/incompressible/windEnergy/ABLSolver/readABLProperties.H

@@ -93,12 +93,9 @@ int main(int argc, char *argv[])
     // field.
     U.correctBoundaryConditions();
     phi = linearInterpolate(U) & mesh.Sf();
+    #include "turbulenceCorrect.H"
     T.correctBoundaryConditions();
   //p_rgh.correctBoundaryConditions();
-    turbulence->correct();
-    Rwall.correctBoundaryConditions();
-    qwall.correctBoundaryConditions();
-
 
     while (runTime.loop())
     {
@@ -113,14 +110,20 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+            Info << "   Predictor..." << endl;
             #include "UEqn.H"
+            #include "turbulenceCorrect.H"
             #include "TEqn.H"
 
             // --- Pressure corrector loop
+            int corr = 0;
             while (pimple.correct())
             {
+                Info << "   Corrector Step " << corr << "..." << endl;
                 #include "pEqn.H"
+                #include "turbulenceCorrect.H"
                 #include "TEqn.H"
+                corr++;
             }
 
             // --- Compute the velocity flux divergence
@@ -130,15 +133,15 @@ int main(int argc, char *argv[])
             #include "correctGradP.H"
 
             // --- Update the turbulence fields
-            if (pimple.turbCorr())
-            {
-                turbulence->correct();
-            }
+//          if (pimple.turbCorr())
+//          {
+//              turbulence->correct();
+//          }
 
             // --- Update the boundary momentum and
             //     temperature flux conditions
-            Rwall.correctBoundaryConditions();
-            qwall.correctBoundaryConditions();
+//          Rwall.correctBoundaryConditions();
+//          qwall.correctBoundaryConditions();
         }   
 
 

applications/solvers/incompressible/windEnergy/ABLTerrainSolver/ABLTerrainSolver.C

@@ -5,7 +5,10 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I../commonAlgorithms \
+    -I../commonAlgorithms/interpolate2D \
+    -I../commonAlgorithms/windRoseToCartesian
 
 
 EXE_LIBS = \

applications/solvers/incompressible/windEnergy/ABLTerrainSolver/Make/options

@@ -2,6 +2,8 @@
 
     #include "computeCoriolisForce.H"
 
+    #include "computeBuoyancyTerm.H"
+
     fvVectorMatrix UEqn
     (
         fvm::ddt(U)                             // time derivative
@@ -24,7 +26,7 @@
             (
                 (
                   - fvc::snGrad(p_rgh)          // modified pressure gradient
-                  - ghf*fvc::snGrad(rhok)       // buoyancy force
+                  + buoyancyTerm                // buoyancy force
                 ) * mesh.magSf()
             )
         );