Some updates to ABLTerrainSolver.

Author: mchurchf

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

@@ -49,7 +49,10 @@
             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())));
 

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

@@ -222,6 +222,29 @@
        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

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

@@ -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/Make/options

@@ -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/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()
             )
         );