Added an option to enable constrainHbyA.

dafoam/pyDAFoam.py

@@ -11,7 +11,7 @@
 
 """
 
-__version__ = "3.0.7"
+__version__ = "3.0.8"
 
 import subprocess
 import os
@@ -491,6 +491,13 @@ def __init__(self):
         ## multiple times to make sure the AD seeds are propagated properly for the iterative BCs.
         self.hasIterativeBC = False
 
+        ## whether to use the constrainHbyA in the pEqn. The DASolvers are similar to OpenFOAM's native
+        ## solvers except that we directly compute the HbyA term without any constraints. In other words,
+        ## we comment out the constrainHbyA line in the pEqn. However, some cases may diverge without
+        ## the constrainHbyA, e.g., the MRF cases with the SST model. Here we have an option to add the
+        ## constrainHbyA back to the primal and adjoint solvers.
+        self.useConstrainHbyA = False
+
         # *********************************************************************************************
         # ************************************ Advance Options ****************************************
         # *********************************************************************************************

src/adjoint/DAResidual/DAResidualRhoSimpleCFoam.C

@@ -142,14 +142,17 @@ void DAResidualRhoSimpleCFoam::calcResiduals(const dictionary& options)
     // copied and modified from pEqn.H
     volScalarField rAU(1.0 / UEqn.A());
     volScalarField rAtU(1.0 / (1.0 / rAU - UEqn.H1()));
-    //volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    volVectorField HbyA(constrainHbyA(rAU * UEqn.H(), U_, p_));
     //***************** NOTE *******************
-    // constrainHbyA has been used since OpenFOAM-v1606; however, We do NOT use the constrainHbyA
-    // function in DAFoam because we found it significantly degrades the accuracy of shape derivatives.
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives.
     // Basically, we should not constrain any variable because it will create discontinuity.
-    // Instead, we use the old implementation used in OpenFOAM-3.0+ and before
-    volVectorField HbyA("HbyA", U_);
-    HbyA = rAU * UEqn.H();
+    // Here we have a option to use the old implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    label useConstrainHbyA = daOption_.getOption<label>("useConstrainHbyA");
+    if (!useConstrainHbyA)
+    {
+        HbyA = rAU * UEqn.H();
+    }
+
     tUEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho_) * fvc::flux(HbyA));
@@ -191,7 +194,6 @@ void DAResidualRhoSimpleCFoam::calcResiduals(const dictionary& options)
     // copied and modified from pEqn.H
     phiRes_ = phiHbyA + pEqn.flux() - phi_;
     normalizePhiResiduals(phiRes);
-
 }
 
 void DAResidualRhoSimpleCFoam::updateIntermediateVariables()

src/adjoint/DAResidual/DAResidualRhoSimpleFoam.C

@@ -177,14 +177,17 @@ void DAResidualRhoSimpleFoam::calcResiduals(const dictionary& options)
     // copied and modified from pEqn.H
     volScalarField rAU(1.0 / UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho_ * rAU));
-    //volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    volVectorField HbyA(constrainHbyA(rAU * UEqn.H(), U_, p_));
     //***************** NOTE *******************
-    // constrainHbyA has been used since OpenFOAM-v1606; however, We do NOT use the constrainHbyA
-    // function in DAFoam because we found it significantly degrades the accuracy of shape derivatives.
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives.
     // Basically, we should not constrain any variable because it will create discontinuity.
-    // Instead, we use the old implementation used in OpenFOAM-3.0+ and before
-    volVectorField HbyA("HbyA", U_);
-    HbyA = rAU * UEqn.H();
+    // Here we have a option to use the old implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    label useConstrainHbyA = daOption_.getOption<label>("useConstrainHbyA");
+    if (!useConstrainHbyA)
+    {
+        HbyA = rAU * UEqn.H();
+    }
+
     tUEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho_) * fvc::flux(HbyA));

src/adjoint/DAResidual/DAResidualSimpleFoam.C

@@ -127,14 +127,17 @@ void DAResidualSimpleFoam::calcResiduals(const dictionary& options)
     scalar pRefValue = 0.0;
 
     volScalarField rAU(1.0 / UEqn.A());
-    //volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U_, p_));
+    volVectorField HbyA(constrainHbyA(rAU * UEqn.H(), U_, p_));
     //***************** NOTE *******************
-    // constrainHbyA has been used since OpenFOAM-v1606; however, We do NOT use the constrainHbyA
-    // function in DAFoam because we found it significantly degrades the accuracy of shape derivatives.
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives.
     // Basically, we should not constrain any variable because it will create discontinuity.
-    // Instead, we use the old implementation used in OpenFOAM-3.0+ and before
-    volVectorField HbyA("HbyA", U_);
-    HbyA = rAU * UEqn.H();
+    // Here we have a option to use the old implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    label useConstrainHbyA = daOption_.getOption<label>("useConstrainHbyA");
+    if (!useConstrainHbyA)
+    {
+        HbyA = rAU * UEqn.H();
+    }
+
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 

src/adjoint/DAResidual/DAResidualSimpleTFoam.C

@@ -146,14 +146,16 @@ void DAResidualSimpleTFoam::calcResiduals(const dictionary& options)
     scalar pRefValue = 0.0;
 
     volScalarField rAU(1.0 / UEqn.A());
-    //volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U_, p_));
+    volVectorField HbyA(constrainHbyA(rAU * UEqn.H(), U_, p_));
     //***************** NOTE *******************
-    // constrainHbyA has been used since OpenFOAM-v1606; however, We do NOT use the constrainHbyA
-    // function in DAFoam because we found it significantly degrades the accuracy of shape derivatives.
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives.
     // Basically, we should not constrain any variable because it will create discontinuity.
-    // Instead, we use the old implementation used in OpenFOAM-3.0+ and before
-    volVectorField HbyA("HbyA", U_);
-    HbyA = rAU * UEqn.H();
+    // Here we have a option to use the old implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    label useConstrainHbyA = daOption_.getOption<label>("useConstrainHbyA");
+    if (!useConstrainHbyA)
+    {
+        HbyA = rAU * UEqn.H();
+    }
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 

src/adjoint/DASolver/DARhoSimpleCFoam/pEqnRhoSimpleC.H

@@ -2,14 +2,17 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0 / UEqn.A());
 volScalarField rAtU(1.0 / (1.0 / rAU - UEqn.H1()));
-//volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+volVectorField HbyA(constrainHbyA(rAU * UEqn.H(), U, p));
 //***************** NOTE *******************
-// constrainHbyA has been used since OpenFOAM-v1606; however, We do NOT use the constrainHbyA
-// function in DAFoam because we found it significantly degrades the accuracy of shape derivatives.
+// constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives.
 // Basically, we should not constrain any variable because it will create discontinuity.
-// Instead, we use the old implementation used in OpenFOAM-3.0+ and before
-volVectorField HbyA("HbyA", U);
-HbyA = rAU * UEqn.H();
+// Here we have a option to use the old implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+label useConstrainHbyA = daOptionPtr_->getOption<label>("useConstrainHbyA");
+if (!useConstrainHbyA)
+{
+    HbyA = rAU * UEqn.H();
+}
+
 tUEqn.clear();
 
 bool closedVolume = false;

src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H

@@ -1,13 +1,16 @@
 volScalarField rAU(1.0 / UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho* rAU));
-//volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+volVectorField HbyA(constrainHbyA(rAU * UEqn.H(), U, p));
 //***************** NOTE *******************
-// constrainHbyA has been used since OpenFOAM-v1606; however, We do NOT use the constrainHbyA
-// function in DAFoam because we found it significantly degrades the accuracy of shape derivatives.
+// constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives.
 // Basically, we should not constrain any variable because it will create discontinuity.
-// Instead, we use the old implementation used in OpenFOAM-3.0+ and before
-volVectorField HbyA("HbyA", U);
-HbyA = rAU * UEqn.H();
+// Here we have a option to use the old implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+label useConstrainHbyA = daOptionPtr_->getOption<label>("useConstrainHbyA");
+if (!useConstrainHbyA)
+{
+    HbyA = rAU * UEqn.H();
+}
+
 tUEqn.clear();
 
 bool closedVolume = false;

src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H

@@ -1,13 +1,15 @@
 {
     volScalarField rAU(1.0 / UEqn.A());
-    //volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    volVectorField HbyA(constrainHbyA(rAU * UEqn.H(), U, p));
     //***************** NOTE *******************
-    // constrainHbyA has been used since OpenFOAM-v1606; however, We do NOT use the constrainHbyA
-    // function in DAFoam because we found it significantly degrades the accuracy of shape derivatives. 
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives.
     // Basically, we should not constrain any variable because it will create discontinuity.
-    // Instead, we use the old implementation used in OpenFOAM-3.0+ and before
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU * UEqn.H();
+    // Here we have a option to use the old implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    label useConstrainHbyA = daOptionPtr_->getOption<label>("useConstrainHbyA");
+    if (!useConstrainHbyA)
+    {
+        HbyA = rAU * UEqn.H();
+    }
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 
     MRF.makeRelative(phiHbyA);

src/adjoint/DASolver/DASimpleTFoam/pEqnSimpleT.H

@@ -1,13 +1,15 @@
 {
     volScalarField rAU(1.0 / UEqn.A());
-    //volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+    volVectorField HbyA(constrainHbyA(rAU * UEqn.H(), U, p));
     //***************** NOTE *******************
-    // constrainHbyA has been used since OpenFOAM-v1606; however, We do NOT use the constrainHbyA
-    // function in DAFoam because we found it significantly degrades the accuracy of shape derivatives. 
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives.
     // Basically, we should not constrain any variable because it will create discontinuity.
-    // Instead, we use the old implementation used in OpenFOAM-3.0+ and before
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU * UEqn.H();
+    // Here we have a option to use the old implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    label useConstrainHbyA = daOptionPtr_->getOption<label>("useConstrainHbyA");
+    if (!useConstrainHbyA)
+    {
+        HbyA = rAU * UEqn.H();
+    }
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 
     MRF.makeRelative(phiHbyA);