Added an option to use constrainHbyA (mdolab#486)

* Added an option to enable constrainHbyA.

* Added the missing header file.

* Used pointer for HbyA.

* Fixed a minor issue.

* Updated the integration test due to an issue for SLSQP v2 scripts.

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/DAResidual.H

@@ -22,6 +22,7 @@
 #include "DAField.H"
 #include "DAFvSource.H"
 #include "IOMRFZoneListDF.H"
+#include "constrainHbyA.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/adjoint/DAResidual/DAResidualRhoSimpleCFoam.C

@@ -142,14 +142,23 @@ 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));
     //***************** 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.
-    // 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();
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives
+    // because constraining variables will create discontinuity. Here we have a option to use the old
+    // implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    autoPtr<volVectorField> HbyAPtr = nullptr;
+    label useConstrainHbyA = daOption_.getOption<label>("useConstrainHbyA");
+    if (useConstrainHbyA)
+    {
+        HbyAPtr.reset(new volVectorField(constrainHbyA(rAU * UEqn.H(), U_, p_)));
+    }
+    else
+    {
+        HbyAPtr.reset(new volVectorField("HbyA", U_));
+        HbyAPtr() = rAU * UEqn.H();
+    }
+    volVectorField& HbyA = HbyAPtr();
+
     tUEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho_) * fvc::flux(HbyA));
@@ -191,7 +200,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,23 @@ 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));
     //***************** 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.
-    // 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();
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives
+    // because constraining variables will create discontinuity. Here we have a option to use the old
+    // implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    autoPtr<volVectorField> HbyAPtr = nullptr;
+    label useConstrainHbyA = daOption_.getOption<label>("useConstrainHbyA");
+    if (useConstrainHbyA)
+    {
+        HbyAPtr.reset(new volVectorField(constrainHbyA(rAU * UEqn.H(), U_, p_)));
+    }
+    else
+    {
+        HbyAPtr.reset(new volVectorField("HbyA", U_));
+        HbyAPtr() = rAU * UEqn.H();
+    }
+    volVectorField& HbyA = HbyAPtr();
+
     tUEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho_) * fvc::flux(HbyA));

src/adjoint/DAResidual/DAResidualSimpleFoam.C

@@ -127,14 +127,22 @@ void DAResidualSimpleFoam::calcResiduals(const dictionary& options)
     scalar pRefValue = 0.0;
 
     volScalarField rAU(1.0 / UEqn.A());
-    //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.
-    // 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();
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives
+    // because constraining variables will create discontinuity. Here we have a option to use the old
+    // implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    autoPtr<volVectorField> HbyAPtr = nullptr;
+    label useConstrainHbyA = daOption_.getOption<label>("useConstrainHbyA");
+    if (useConstrainHbyA)
+    {
+        HbyAPtr.reset(new volVectorField(constrainHbyA(rAU * UEqn.H(), U_, p_)));
+    }
+    else
+    {
+        HbyAPtr.reset(new volVectorField("HbyA", U_));
+        HbyAPtr() = rAU * UEqn.H();
+    }
+    volVectorField& HbyA = HbyAPtr();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 

src/adjoint/DAResidual/DAResidualSimpleTFoam.C

@@ -146,14 +146,22 @@ void DAResidualSimpleTFoam::calcResiduals(const dictionary& options)
     scalar pRefValue = 0.0;
 
     volScalarField rAU(1.0 / UEqn.A());
-    //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.
-    // 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();
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives
+    // because constraining variables will create discontinuity. Here we have a option to use the old
+    // implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    autoPtr<volVectorField> HbyAPtr = nullptr;
+    label useConstrainHbyA = daOption_.getOption<label>("useConstrainHbyA");
+    if (useConstrainHbyA)
+    {
+        HbyAPtr.reset(new volVectorField(constrainHbyA(rAU * UEqn.H(), U_, p_)));
+    }
+    else
+    {
+        HbyAPtr.reset(new volVectorField("HbyA", U_));
+        HbyAPtr() = rAU * UEqn.H();
+    }
+    volVectorField& HbyA = HbyAPtr();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 

src/adjoint/DASolver/DARhoSimpleCFoam/pEqnRhoSimpleC.H

@@ -2,14 +2,23 @@ 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));
 //***************** 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.
-// 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();
+// constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives
+// because constraining variables will create discontinuity. Here we have a option to use the old 
+// implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+autoPtr<volVectorField> HbyAPtr = nullptr;
+label useConstrainHbyA = daOptionPtr_->getOption<label>("useConstrainHbyA");
+if (useConstrainHbyA)
+{
+    HbyAPtr.reset(new volVectorField(constrainHbyA(rAU * UEqn.H(), U, p)));
+}
+else
+{
+    HbyAPtr.reset(new volVectorField("HbyA", U));
+    HbyAPtr() = rAU * UEqn.H();
+}
+volVectorField& HbyA = HbyAPtr();
+
 tUEqn.clear();
 
 bool closedVolume = false;

src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H

@@ -1,13 +1,22 @@
 volScalarField rAU(1.0 / UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho* rAU));
-//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.
-// 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();
+// constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives
+// because constraining variables will create discontinuity. Here we have a option to use the old 
+// implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+autoPtr<volVectorField> HbyAPtr = nullptr;
+label useConstrainHbyA = daOptionPtr_->getOption<label>("useConstrainHbyA");
+if (useConstrainHbyA)
+{
+    HbyAPtr.reset(new volVectorField(constrainHbyA(rAU * UEqn.H(), U, p)));
+}
+else
+{
+    HbyAPtr.reset(new volVectorField("HbyA", U));
+    HbyAPtr() = rAU * UEqn.H();
+}
+volVectorField& HbyA = HbyAPtr();
+
 tUEqn.clear();
 
 bool closedVolume = false;

src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H

@@ -1,13 +1,21 @@
 {
     volScalarField rAU(1.0 / UEqn.A());
-    //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. 
-    // 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();
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives
+    // because constraining variables will create discontinuity. Here we have a option to use the old
+    // implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    autoPtr<volVectorField> HbyAPtr = nullptr;
+    label useConstrainHbyA = daOptionPtr_->getOption<label>("useConstrainHbyA");
+    if (useConstrainHbyA)
+    {
+        HbyAPtr.reset(new volVectorField(constrainHbyA(rAU * UEqn.H(), U, p)));
+    }
+    else
+    {
+        HbyAPtr.reset(new volVectorField("HbyA", U));
+        HbyAPtr() = rAU * UEqn.H();
+    }
+    volVectorField& HbyA = HbyAPtr();
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 
     MRF.makeRelative(phiHbyA);

src/adjoint/DASolver/DASimpleTFoam/pEqnSimpleT.H

@@ -1,13 +1,21 @@
 {
     volScalarField rAU(1.0 / UEqn.A());
-    //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. 
-    // 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();
+    // constrainHbyA has been used since OpenFOAM-v1606; however, it may degrade the accuracy of derivatives
+    // because constraining variables will create discontinuity. Here we have a option to use the old
+    // implementation in OpenFOAM-3.0+ and before (no constraint for HbyA)
+    autoPtr<volVectorField> HbyAPtr = nullptr;
+    label useConstrainHbyA = daOptionPtr_->getOption<label>("useConstrainHbyA");
+    if (useConstrainHbyA)
+    {
+        HbyAPtr.reset(new volVectorField(constrainHbyA(rAU * UEqn.H(), U, p)));
+    }
+    else
+    {
+        HbyAPtr.reset(new volVectorField("HbyA", U));
+        HbyAPtr() = rAU * UEqn.H();
+    }
+    volVectorField& HbyA = HbyAPtr();
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 
     MRF.makeRelative(phiHbyA);

tests/refs/DAFoam_Test_IntegrationRef.txt

@@ -1,23 +1,23 @@
 Dictionary Key: alpha
-@value    4.983717286274688 1e-06 1e-08
+@value    4.968905604685552 1e-06 1e-08
 Dictionary Key: shapey
-@value -0.007921636479688784 1e-06 1e-08
-@value -0.00792163647968879 1e-06 1e-08
-@value 0.007921636479688784 1e-06 1e-08
-@value 0.007921636479688784 1e-06 1e-08
-@value -0.03828250734496957 1e-06 1e-08
-@value -0.03828250734496957 1e-06 1e-08
-@value  0.01296166604998152 1e-06 1e-08
-@value  0.01296166604998153 1e-06 1e-08
-@value -0.01007603172079106 1e-06 1e-08
-@value -0.01007603172079106 1e-06 1e-08
-@value 0.006776111843578327 1e-06 1e-08
-@value 0.006776111843578336 1e-06 1e-08
-@value -0.006275224171581452 1e-06 1e-08
-@value -0.00627522417158145 1e-06 1e-08
-@value 0.006469681351950167 1e-06 1e-08
-@value  0.00646968135195016 1e-06 1e-08
-@value -0.001009411824693773 1e-06 1e-08
-@value -0.001009411824693775 1e-06 1e-08
-@value 0.001009411824693774 1e-06 1e-08
-@value  0.00100941182469377 1e-06 1e-08
+@value 8.655069785387594e-05 1e-06 1e-08
+@value 8.655069785387594e-05 1e-06 1e-08
+@value -8.655069785387605e-05 1e-06 1e-08
+@value -8.655069785387605e-05 1e-06 1e-08
+@value -0.01119217936982357 1e-06 1e-08
+@value -0.01119217936982357 1e-06 1e-08
+@value -0.007799868290607725 1e-06 1e-08
+@value -0.007799868290607725 1e-06 1e-08
+@value 0.0004216141178557599 1e-06 1e-08
+@value 0.0004216141178557599 1e-06 1e-08
+@value -0.001760513948823865 1e-06 1e-08
+@value -0.001760513948823865 1e-06 1e-08
+@value 0.0005508325823005192 1e-06 1e-08
+@value 0.0005508325823005192 1e-06 1e-08
+@value 0.001951434242102298 1e-06 1e-08
+@value 0.001951434242102298 1e-06 1e-08
+@value 9.201202606960528e-05 1e-06 1e-08
+@value 9.201202606960528e-05 1e-06 1e-08
+@value -9.201202606960528e-05 1e-06 1e-08
+@value -9.201202606960528e-05 1e-06 1e-08

tests/runTests_Integration.py

@@ -104,9 +104,9 @@
 }
 
 optOptions = {
-    "ACC": 1.0e-5,  # convergence accuracy
-    "MAXIT": 2,  # max optimization iterations
-    "IFILE": "opt_SLSQP.out",
+    "tol": 1.0e-5,
+    "max_iter": 2,
+    "output_file": "opt_IPOPT.txt",
 }
 
 # DVGeo
@@ -217,8 +217,8 @@ def setMultiPointObjFuncsSens(xDVs, funcsMP, funcsSens, funcsSensMP, index):
 if gcomm.rank == 0:
     print(optProb)
 
-opt = OPT("slsqp", options=optOptions)
-histFile = os.path.join("./", "slsqp_hist.hst")
+opt = OPT("ipopt", options=optOptions)
+histFile = os.path.join("./", "ipopt_hist.hst")
 sol = opt(optProb, sens=optFuncs.calcObjFuncSensMP, storeHistory=histFile)
 
 if gcomm.rank == 0: